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INTRODUCTION 


This  report  presents  a  study  of  the  indirect  boundary  eleaent 
■ethod  and  its  potential  advantages  for  solving  one-  and  two- dimensional 
linear  structural/stress  analysis  problems.  Currently,  the  finite 
eleaent  method  is  very  capable  of  performing  these  tasks.  However,  the 
boundary  eleaent  method  potentially  offers  an  opportunity  for  increased 
productivity  in  these  areas. 

Background 

Because  the  application  of  the  boundary  eleaent  method  requires 
only  that  the  boundary  of  the  structure  be  subdivided,  as  contrasted 
with  the  requirement  that  the  entire  domain  of  the  structure  be  sub¬ 
divided  when  applying  the  finite  eleaent  method,  the  boundary  element 
method  may  increase  productivity  in  linear  structural  analysis.  Some¬ 
times  the  former  method  is  called  a  boundary  method  and  the  latter  a 
domain  method. 

A  good  account  of  boundary  eleaent  methods  from  the  perspective  of 
finite  eleaent  methods  can  be  found  in  Reference  1.  There  are  two 
types  of  boundary  eleaent  methods:  a  direct  method  and  an  indirect 
method.  The  difference  between  the  two  methods  is  not  easily  explained 
in  a  brief  manner.  The  indirect  method  is  rather  intuitive,  while  the 
direct  method  is  more  formal.  Both  methods  are  equally  effective  in 
general.  A  comparison  of  indirect  and  direct  boundary  eleaent  methods 
can  be  found  in  Reference  2.  An  early  treatment  of  indirect  boundary 
methods  can  be  found  in  Reference  3. 

Accuracy  aside  for  the  moment,  the  effectiveness  of  any  numerical 
stress  analysis  procedure  depends  primarily  on  the  manual  effort  required 
for  pre-  and  post-processing  of  the  required  input  and  output  data,  and 
to  a  lesser  extent  on  the  computer  usage  cost  of  the  associated  structural 


computer  program.  The  boundary  element  method  requires  considerably 
less  input  and  output  data  preparation,  particularly  input  data,  because 
fewer  subdivisions  are  necessary  to  describe  the  structure  boundary  than 
the  structure  itself.  The  manual  effort  associated  with  pre-  and  post¬ 
processing  finite  element  data  is  often  very  considerable  and  frus- 
tratingly  long,  notwithstanding  the  advantages  of  current  automated 
techniques.  (Obviously,  such  techniques  are  also  applicable  to  the 
boundary  element  methods  as  well.)  Thus,  the  boundary  element  methods 
offer  potential  savings  by  reducing  the  manual  labor  of  the  stress 
analyst  and  the  engineering  technician,  particularly  at  the  input  data 
level,  for  a  given  structural  problem. 

The  second  advantage  of  the  boundary  element  method  involves  reduced 
computational  effort  in  the  structural  analysis  computer  program.  In 
linear,  static  structural  analysis  by  the  finite  element  method  it  is 
well  known  that  most  of  the  computational  cost  lies  in  the  solution  of 
the  system  of  linear  algebraic  equations  that  result  from  the  finite 
element  subdivision.  In  finite  element  computer  programs  very  effective 
Gaussian  elimination,  equation-solving  algorithms  have  evolved  that 
minimize  this  cost  (Ref  4).  However,  for  the  same  structural  analysis 
problem,  the  number  of  linear  algebraic  equations  that  must  be  solved 
when  employing  the  boundary  element  method  is  generally  far  smaller  than 
in  the  finite  element  method.  This  is  because  the  boundary  element 
approach  immediately  reduces  the  structural  problem  by  one  dimension  due 
to  the  necessity  of  having  only  to  subdivide  the  boundary  of  the  struc¬ 
ture.  Thus,  three-dimensional  problems,  as  in  the  stress  analysis  of 
solids,  are  reduced  to  two-dimensional  problems;  two-dimensional  problems, 
as  in  the  stress  analysis  of  membranelike  plates,  are  reduced  to  one-dimensional 
problems;  and  one-dimensional  problems,  as  in  the  analysis  of  beams,  are 
reduced  to  what  can  be  termed  "point  problems."  In  theory,  this  amounts 
to  a  distinct  computational  advantage  of  the  boundary  element  method 
over  the  ubiquitous  finite  element  method.  There  are  mitigating 
considerations,  however.  In  the  case  of  a  materially  homogeneous  problem, 
for  example,  the  coefficient  matrix  in  the  linear  algebraic  system  is 
full  in  the  boundary  eleittnt  method,  whereas  the  coefficient  matrix, 
though  much  larger,  is  both  sparse  and  symmetric  in  the  finite  element 
method. 
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The  objective  of  this  study  is  to  assess  the  accuracy  and  potential 
of  the  indirect  boundary  element  method  in  linear  structural  analysis 
through  numerical  experimentation.  The  indirect  boundary  element  method 
is  to  be  explicated  and  then  demonstrated  by  developing  a  one-dimensional 
computer  program  and  a  two-dimensional  computer  program.  Another  objec¬ 
tive  is  to  determine  the  suitability  of  the  method  as  a  structural/stress 
analysis  tool  when  implemented  on  microcomputers. 

Scope 


The  theoretical  formulation  of  the  indirect  boundary  element  method 
is  illustrated  first  by  developing  the  framework  of  one-dimensional 
beams  resting  on  elastic  foundations,  and  then  extending  the  same  concept 
to  the  framework  of  two-dimensional  plane  stress  or  plane  strain  elasto- 
statics. 

Computer  programs  are  written  both  in  BASIC  and  FORTRAN  that  numer¬ 
ically  implement  the  theoretical  formulations  for  the  one-dimensional 
application.  The  program  for  the  two-dimensional  application  is  written 
in  FORTRAN.  Accuracy  of  the  indirect  boundary  element  solution  is 
assessed  through  comparison  with  theoretical  solutions  and  with  solu¬ 
tions  from  the  alternative,  direct  boundary  element  method. 


THEORETICAL  BASIS  OF  THE  INDIRECT  BOUNDARY  ELEMENT  METHOD 

The  indirect  boundary  element  method  is  a  general  numerical  solu¬ 
tion  technique  for  solving  boundary  value  problems  in  engineering  science. 

A  bibliography  is  included  in  this  report  which  shows  the  breadth  of 
engineering  boundary  value  problems  that  can  be  approached  with  the 
method.*  The  necessary  theoretical  relationships  and  equations  for  the 
method  as  applied  to  structural  problems  com  from  the  theory  of  elasticity. 


Also  see  Reference  5  for  a  wide  selection  of  applications. 


This  theoretical  basis  is  explicated  herein,  first  with  a  one-dimensional 
example  of  a  beam  resting  on  an  elastic  foundation,  and  second,  with  the 
general  problem  of  two-dimensional  elastostatics. 

Once  the  boundary  value  problem  has  been  completely  stated,  the 
numerical  solution  of  that  problem  by  the  indirect  boundary  element 
method  follows  three  basic  steps. 

1.  Establish  the  infinite-domain  Green's  functions  appropriate  to 
the  boundary  value  problem. 

2.  Form  and  solve  the  auxiliary  boundary  value  problem  in  the 
infinite  doaiain  by  employing  superposition  of  the  established  solutions. 

3.  Invoke  the  Kirchhoff  uniqueness  theorem  to  obtain  the  solution 
to  the  original  boundary  value  problem  from  the  solution  of  the  aux¬ 
iliary  problem. 

Since  these  steps  also  contain  information  on  the  natural  limitations  of 
the  indirect  boundary  element  method,  they  are  discussed  below. 

A  Green's  function  is  a  known  solution  to  the  governing  differ¬ 
ential  equation  of  the  given  boundary  value  problem.  It  is  very  much 
like  an  influence  function,  a  concept  familiar  to  undergraduate  civil 
engineers,  which  algebraically  determines  the  response,  say  displacement 
at  some  field  point,  due  to  a  prescribed  unit  concentrated  force  at  some 
other  point,  called  a  source  point.  The  important  concept  here  is  that 
such  a  solution  to  the  governing  differential  equation  must  be  known  at 
the  outset  for  the  method  to  be  applicable.*  Step  1  implies  that  the 
appropriate  Green's  functions  must  exist. 

The  key  word  in  Step  2  is  superposition.  Thus,  the  boundary  value 
problem  must  be  linear  for  the  indirect  boundary  element  method  to  be 
applicable  as  a  solution  technique.  It  should  be  noted,  however,  that 
despite  this  limitation,  some  nonlinear  problems  are  solved  with  boundary 


It  is  also  true  that  a  Green's  function  satisfies  boundary  conditions 
as  well.  See  References  6  and  7  for  good  accounts  of  Green's  functions. 
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element  Methods  (Ref  5  and  8) .  These  approaches  must  inevitably  use  a 
sequence  of  linearizations.  Nonetheless,  the  principal  application  of 
the  boundary  eleaent  Method  at  present  is  to  linear  problems.  The 
solution  to  the  auxiliary  problem  is  built  up  from  the  superposition  of 
unit  solution  coaponents  provided  by  the  known  Green's  functions. 

Finally,  Step  3  implies  that  the  solution  to  the  actual  problem  can 
be  obtained  only  if  Kirchhoff's  uniqueness  theorem  can  be  invoked. 
According  to  Reference  9,  this  theorem  states: 

If,  in  addition  to  the  body  forces,  either  the  surface 
forces  or  the  surface  displacements  are  given  on  the  boundary 
of  an  elastic  body,  there  exists  only  one  form  of  equilibrium 
in  the  sense  that  the  distribution  of  stresses  and  strains  in 
the  body  is  determined  uniquely. 

The  theorea  requires  that  the  structural  problea  be  limited  to  infinites 
iaal  strains  and  displaceswnts .  Exactly  why  this  theorem  is  invoked 
will  becoae  clearer  when  applications  of  the  indirect  boundary  element 
method  are  presented  in  the  next  section. 

Beam  Resting  on  an  Elastic  Foundation 

The  necessary  equations  to  be  programmed  into  a  computer  for  the 
nuaerical  solution  of  beams  on  elastic  foundations  by  the  indirect 
boundary  eleaent  are  developed  below.  This  development  closely  follows 
the  account  given  in  Reference  10.  The  class  of  problems  addressed  is 
shown  in  Figure  1.  The  boundary  value  problem  is  stated  mathematically 
as  follows.  Solve  the  fourth  order  differential  equation 

.4 

El  =  b(x)  -  k  u(x)  (1 

dx* 

where  u  =  lateral  deflection 


b  =  prescribed  lateral  load 

k  =  elastic  foundation  stiffness  per  unit  length 
El  =  beaa  bendi,  —’gid  * 


The  equation  is  to  be  solved  subject  to  the  following  four  boundary 
conditions : 


at  x  = 

o, 

o 

ii 

8 

and 

s  =  0 

at  x  = 

L, 

0  =  0 

and 

s  =  0 

where  m  =  bending  moment 
s  =  shear  force 
8  =  slope 


(2) 


The  boundary  condition  at  x  =  0  implies  a  free  end  condition,  and 
at  x  =  L  a  symmetry  condition  is  implied.  Other  beam  boundary  conditions 
can  also  be  imposed. 

Step  1  of  the  procedure  requires  the  Green's  functions  to  be  given 
for  the  above  problem.  These  functions  can  be  found  in  Reference  11  and 
are  applicable  to  an  infinitely  long  beam  and  foundation  as  shown  in 
Figure  2.  The  appropriate  beam  response  functions  at  a  field  point  Q 
for  a  unit  concentrated  force  at  point  P  are 


i(r) 

=  Ae 

Pr  (cos  0r  +  sin  0r) 
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=  e~ 
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where  sgn(y-x)  =  \ 

-1  y  <  x 

1 

undefined  y  =  x 

and 

P  = 

Sjk/U El 

(3) 


Similarly,  the  required  Green's  functions  for  the  beam  response  at  Q  due 
to  a  unit  moment  at  P  are 


u  (r)  = 


0(r)  =  e  sin  pr  •  sgn(y-x) 


6*(r) 


-  P  m(r)  =  -  J—  e  (cos  {Jr  -  sin  Pr) 


m*(r) 


6  1 

=  s(r)  =  — —  cos  Pr  •  sgn(y-x) 


s  (r) 


=  -g-u(r)  =  P  ^2 —  (cos  +  sin  ^ 


In  Step  2  the  auxiliary  problem  is  formed  first  by  imbedding  the 
actual  beam  and  foundation  complete  with  loading  into  the  infinite 
domain  as  shown  in  Figure  3.  It  is  apparent  that  the  required  boundary 
conditions  at  points  1  and  2  are  not  imposed  in  the  infinite  beam  unless 
something  else  is  done  to  enforce  them.  Therefore,  the  unknown  forces  iji 
and  unknown  moments  tpt,  are  applied  to  these  two  points  to  impose  the 
given  boundary  conditions.  The  auxiliary  problem  can  now  be  stated  as: 
find  the  unknown  forces  and  moments  at  points  1  and  2  such  that  the 
prescribed  boundary  conditions  of  the  actual  problem  are  satisfied. 

The  distributed  load  b(x)  can  be  resolved  into  N  statically  equiva¬ 
lent  concentrated  forces  acting  on  the  beam,  and,  in  general,  there  may 
also  be  M  concentrated  moments  acting  on  the  actual  beam.  Using  the 
Green's  functions  and  the  principle  of  superposition,  a  set  of  linear 
algebraic  equations  in  the  two  unknown  forces  (|k  and  the  two  unknown 
moments  t|#  can  be  established  such  that  the  required  boundary  conditions 
are  satisfied  upon  their  solution.  Construction  of  this  system  of 
equations  proceeds  as  follows: 


For  Sj  =  0  at  x  =  +e 

*  £  *  -fr  £ 

Sj  =  stejO)^  +  s  (e.OMj  +  s(e,L)«|»2  +  s  (e,L)i|i2 


+  jH  s(e,x,)b.  +  £  s  (e,x.)b.  =  0 


B1  = 


iCe.O)^  +  a  (e.O)^  +  m(e,L)»|>2  +  a  (e,L)»|f2 


N  - 


M 


+  Y  ■(£,*. )b.  +  Y  ®  (e»x4)bi  = 
i=l  1  1  j=i  J  J 


For  s2  =  0  at  x  =  L-e 


*  -k  it 

siL-£. +  s  (L-e +  s(L-e,L)»|»2  +  s  (L-e,L)«|»2 


N  * 


+  Y  s(L-e,x.)b.  +  Y  8  (b-e.xJb. 
i=l  1  1  j=l  J  J 


=  0 


For  ©2  =  0  at  x  =  L-e 


02  = 


ed-e^Hj  +  0*(L-e,O)iJ»*  +  0(L-e,L)«|»2  +  0*(L-£,L)*2 

N  a  M  .  . 

+  Y  0(L-e,x.)b.  +  Y  6  (L-e,x.)b.  =  0 
i=l  1  1  j=l  J  J 


The  ter*  £  represents  an  arbitrarily  small  distance  to  indicate  that 
these  functions  are  evaluated  just  inside  the  actual  or  real  beam  domain 
In  matrix  form,  these  equations  are  written  as 


«(e,0) 

•*(e,0) 

A  * 

«(e,L)  s  (e,L) 

*1 

■(e,0) 

■* (e,0) 

■(e,L)  ■*(£,!) 

' 

★ 

*1 

*(L-e,0) 

•*«.-£, 0) 

»(L-c,L) 

* 2 

A 

.e(L-e,0) 

6*(L-e,0) 

e(L-c,L)  e*(l-£,L) 

< 

•(«.*!) 

•(c,*2) 

A  1  * 

...  •(fi.Xj,)  |  t  (£, 

■u,x2> 

...  |  ■*(£, 

*l) 

«a-c,Xj) 

•(L-c,x2) 

...  »(L-£,Xj|)  |  »*(L-£,Xj) 

.ea-c,x  ) 

6(L-e,x2) 

...  ea-S'X,.)  I  e*(L-E,x.) 

♦ 


•*(e,x2) 

**(L-e,x2) 

6*(L-c,x,) 


(5) 

T 


0 


or 


G  *  ♦  H  b  =  F 


(5a) 


The  solution  to  the  auxiliary  problem  is  then  given  by 
*  =  G_1 (F  -  H  b) 

At  M  'a#  m  a/' 


(6) 


m 


& 


§ 


It  is  interesting  to  note  that  the  order  of  the  coefficient  matrix 
G  is  only  four,  and  that  this  will  always  be  the  case  independently  of 
the  beam  and  foundation  length  so  long  as  the  they  are  prescribed  homo* 
geneous.  Thus,  the  solution  will  always  require  that  matrices  of  only 
order  four  be  inverted.  The  one-dimensional  boundary  value  problem  has 
been  reduced  to  a  discrete  boundary  element  problem  involving  only  two 
points  in  the  domain  -  the  boundary  points  1  and  2  of  the  actual  beam 
and  foundation  problem.  This  is  in  contrast  to  a  standard  finite  element 
solution  approach  that  would  have  reduced  the  continuous  problem  to  a 
discrete  problem  at  n  points  (nodes)  in  the  domain  along  the  beam. 

For  all  practical  purposes  the  greater  part  of  the  solution  effort 
in  the  indirect  boundary  element  method  is  accomplished  in  Step  2. 

Since  the  forces  and  moments  at  the  end  points  are  now  known  for  the 
auxiliary  problem,  and  since  they,  by  definition,  impose  the  prescribed 
boundary  conditions  in  the  infinite  domain,  Step  3  can  be  addressed. 
Kirchhoff's  uniqueness  theorem  simply  specifies  that  the  solution  in  the 
domain  [0,L]  of  the  actual  problem  is  the  same  as  the  solution  for  the 
dosuin  [0,L]  of  the  auxiliary  problem.  The  solution  of  the  actual 
problem  can  then  be  obtained  by  superposition  of  the  Green's  functions 
applied  to  the  auxiliary  problem  as  follows. 

Suppose  the  beam  response  is  desired  at  some  set  of  K  points  within 
the  dosuin  [0 ,L] .  Then  for  k  =  1,2,...,  K  the  displacesient,  moment, 
slope,  and  shear  force  can  be  written  directly  as 


u(xfc)  =  u(xk,0)ij»1  +  u(xk,L)<t»2  +  u*(xk,0)«j»*  +  u*(xR,L)4»2 


+  £  u(xk,x.)b.  +  u*(x.,x.)b* 

i=l  *  1  1  j=l  K  J  J 


mm 
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*  *. 


»(xk)  =  ■(xk,0)4»1  +  m(xk,L)4»2  +  m  (x^O)^  +  m  (xk,L)t|>2 
N  -  M  *  * 

♦  £  ■<vxi)bi  +  E  * 

1=1  J=1  J  J 


e(xk)  =  ecxj^.oHj  +  e(xk,D«|<2  +  e  (xk,o)«j>j  +  e  (xk,u»i»2 

N  a  H 

*  e(V*i)bi  *  ?  6*(xk,nj)bj 

1=1  j=l  J  J 

A  A  ^  ^  ^ 
s(xk)  =  s(xk,0)«|»1  +  s(xk,L)4»2  +  s  (Xj^O)^  +  s  (xk,U»|»2 

N  a  M 

+  £  s(xk»xi^bi  +  2Z  s*(xk*xi)bt 

i=l  *  1  1  j=l  K  3  3 

In  matrix  fora  these  equations  are 


u(xk) 

’A  *  it 

u(xk,0)  u(xk,I)  u  (xk,0)  u  (*k,t) 

*1 

.(xk) 

■(*k»0)  a(xk,L)  a*(xk,0)  a*(xk,L) 

*2 

*  "  *  * 

★ 

®C*k) 

e(xk,o)  e(xk,L)  e  (xk,o)  e  (xk,u 

*1 

«(xk) 

* 

«(*k,0)  •(xk,t)  •*(*k>0)  •*Cxk,L) 

* 

+2 

u(xk,Xj)  u(*k,x2)  ...  u(xk,Xg)  |  u*(xk,Xj)  u*(xk,x2)  .. 
■(xk,Xj)  ■(xk,x2)  ...  ■(xk,xM)  |  ■*(xk,x1)  a*(xk,x2)  .. 
»Cxk,x1)  e(xk,x2)  ...  e(xk,xlf)  I  e  (xk,xj)  e  (xk,x2)  .. 
•  (xk>Xj)  *(xk>x2)  * •  •  I  ®  ®  ^xk,x2^ 


u*(xk,xM) 

■*  tvV 

e*(*k,*M) 


b* 

b* 


Note  that  no  further  inversion  of  matrices  is  necessary.  The  above 
equations  are  merely  algebraic  equations  giving  the  desired  solution  to 
the  actual  beam  and  foundation  problem  at  a  set  of  prescribed  points. 
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The  rectangular  matrix  is  the  largest  matrix  to  be  formed  and  it  is 
order  4x(N+M),  where  N  is  the  number  of  concentrated  or  statically 
equivalent  lateral  forces  prescribed  on  the  beam,  and  M  is  the  number  of 
concentrated  moments  prescribed  along  the  length  of  the  beam. 

With  the  largest  matrix  to  be  inverted  (see  Step  2)  an  order  four 
matrix,  and  the  largest  matrix  to  be  formed  (but  not  necessarily  stored) 
a  4x(N+M)  matrix,  the  indirect  boundary  method  appears  to  be  ideally 
suited  to  microcomputers. 

TWo-Dimensional  Elastostatics 


The  application  of  the  indirect  boundary  element  method  presented 
here  is  for  the  calculation  of  stresses  and  displacements  in  the  plane 
of  two-dimensional  elastic  plates  that  are  subject  to  stretching  due  to 
prescribed  edge  forces  and  edge  displacements.  The  form  of  this  presen¬ 
tation  follows  that  given  in  Reference  12.  Figure  4  illustrates  the 
general  elastostatics  problem  under  consideration.  The  normal  stress 
components  are  and  ay,  and  the  shear  stress  component  is  x  ,  rela¬ 
tive  to  the  x-y  coordinate  system  shown.  The  corresponding  displacement 
components  are  u  and  v.  The  mechanical  properties  of  the  plate  are 
defined  by  the  aradulus  of  elasticity,  E,  and  Poisson's  ratio,  p. 

The  fonal  boundary  value  problem  can  be  stated  as  the  requirement 
for  a  solution  to  the  two-dimensional  stress  equilibrium  partial  dif¬ 
ferential  equations  (ignoring  body  forces)  in  the  domain  Q: 


3x 

_jar 

8y 


3x_7 

_ 52 

8x 


0 


in  n 


0 


(8) 


And  further,  that  the  solution  also  satisfy  the  stress  and  displacement 
boundary  conditions 


where  Fq  and  F  are  prescribed  tractions,  in  the  normal  and  tangential 
directions,  respectively,  along  a  section  of  the  boundary  Tp,  and  u  and 
v  are  the  prescribed  components  of  displacement  in  the  x  and  y  direc¬ 
tions,  respectively,  along  the  reauining  section  of  the  boundary  rQ. 

The  indirect  boundary  element  approach  to  the  solution  begins  with 
the  establishawnt  of  the  Green's  functions  appropriate  to  the  governing 
partial  differential  equation.  In  this  case,  these  functions  are  solu¬ 
tions  for  the  three  stress  components  and  the  two  displacement  compo¬ 
nents  at  a  field  point  Q  due  to  a  unit  concentrated  force  at  a  source 
point  P  in  the  plane  of  an  infinite  two-dimensional  elastic  region  as 
shown  in  Figure  5.  The  following  functions  can  be  found  in  Reference  12 
and  pertain  to  plane  strain  conditions  (i.e.,  the  strain  component 
nonaal  to  the  flat  plate  is  zero,  and  the  stress  component  normal  to  the 
plate  is,  in  general,  nonzero).  These  functions  are  generally  given  in 
terms  of  polar  coordinates  in  standard  references  on  the  theory  of 
elasticity;  see,  for  example,  Reference  13.  Nonetheless,  the  Green's 
functions  in  the  Cartesian  xy  system  are 


ax(Q,P)  = 


cos  0 


2  n 


2X  +  3G  .  2  (A  +  G)  y2] 
A  +  2G  2 


A  +  2G 


x 

T 


V,’P)  ■ 


u(Q,P)  = 


cos  6 

G  . 

2(A  G)  x2l  v 

2  n 

A  ♦  2G 

A  +  2G  22 

r  J  r 

sin  8 

G  t 

2 (A  +  G)  y2]  x 

2  n 

A  ♦  2G 

A  +  2G  22 

r  J  r 

cos  8  r 

A  +  3G 

In  r  X  +  G 

2  71  L  2G(A  +  2G) 

2G(A  +  2G) 

.  sin  8  f 

A  +  G 

xjlI 

2  n  [2G(A  +  2G) 

r2  J 

sin  6  [ 

A  +  3G 

In  r  A  +  G 

2  7C  L  2G(A  ♦  2G) 

lD  r  2G(A  ♦  2G) 

.  cos  8  r 

A  +  G 

*2.1 

2  71  L2G(A  +  2G) 

7J 

where 


r  =  V7T7 


A  = 


E  U 


(1  +  m)(1  -  2jj) 


G  = 


2(1  +  M) 


(10c) 


(11a) 


(lib) 


(12a) 

(12b) 


(12c) 


Though  these  equations  are  for  plane  strain  conditions,  they  can  easily 
be  converted  to  apply  to  plane  stress  conditions  with  the  substitution 
everywhere  of  A'  for  A  where 


A’ 


2  A  G 
A  +  2G 


(13) 


In  that  event  the  class  of  problems  being  considered  would  be  such  that 
the  stress  component  normal  to  the  plane  is  zero,  and  the  strain  compo¬ 
nent  normal  to  the  plane  is  nonzero. 

Step  2  of  the  indirect  boundary  element  procedure  begins  with  the 
formation  of  the  auxiliary  problem  by  scribing  the  actual  problem  boun¬ 
dary  on  the  infinite  sheet  of  the  same  material  and  thickness.  An 
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unknown  system  of  normal  and  tangential  tractions,  Pq  and  Pt,  are  applied 
to  the  scribed  boundary  in  the  plane.  These  forces  are  shown  in  Figure  6. 
The  specific  task  is  to  solve  the  auxiliary  problem  for  this  unknown 
system  such  that  the  given  boundary  conditions  are  satisfied  in  the 
infinite  sheet. 

The  Green's  functions  (Equations  10  and  11)  and  superposition  are 
employed  to  write  equations  relating  the  stresses  and  displacements  at  a 
field  point  Q  on  the  scribed  boundary  (but  also  in  the  infinite  sheet) 
in  terms  of  the  unknown  forces  at  a  source  point  P  as  P  is  successively 
aioved  around  the  entire  scribed  boundary.  That  is,  the  influence  of  the 
tractions  at  each  source  point  P  on  the  boundary  is  superimposed  at  the 
field  point  Q  in  forming  the  equations  for  the  total  response  at  field 
point  Q.  Since  the  unknown  force  system  is  continuous,  the  superposition 
is  accomplished  by  an  integral  on  the  boundary,  and  a  boundary  integral 
equation  corresponding  to  the  field  point  Q  results. 

Since  the  integration  is  around  the  entire  boundary  £,  there  will 
occur  a  singularity  condition  when  the  source  point  P  and  the  field 
point  Q  coincide.  To  provide  for  this  instance,  the  boundary  integral 
is  divided  into  two  integrals,  one  along  a  small  segment  M  that  contains 
the  field  point  Q,  and  the  other  along  the  remaining  portion  of  the 
boundary  i  -  A£.  The  contribution  of  the  former  integral  is  evaluated 
in  the  limit  as  the  distance  r  between  P  and  Q  vanishes.  Determining 
the  singularity  value  of  the  integrals  is  a  lengthy  mathematical  exercise 
that  is  omitted  here.  However,  this  singularity  contribution  is  important 
as  it  will  form  block  diagonal  entries  (or  submatrices)  in  the  coefficient 
matrix  of  the  system  of  algebriac  equations  that  will  result  from  the 
numerical  implementation . 

The  boundary  integral  equations  thusly  formed  for  an  arbitrary 
boundary  point  Q  are  given  below  in  Equations  14  and  15,  where  the  first 
terms  on  the  right-hand  side  in  brackets  are  contributed  by  the 
singularity  condition. 
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where  c 
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1  + 


2(A  +  G) 
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G _ 

C1  "  2n(A  +  2G) 

_  A  +  G 

C2  n(A  +  2G) 

A  +  3G 

c3  4nG(A  +  2G) 
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At  an  arbitrary  point  Q  on  the  boundary,  the  prescribed  boundary 
condition  either  for  stress  or  for  displacement  is  imposed.  The  unknown 
force  distributions  PQ  and  Pt  appearing  on  the  right-hand  sides  of 
Equations  14  or  15  will  correspond  to  these  conditions.  The  equilibrium 
of  a  material  point  Q  on  Tp  is  considered  with  the  aid  of  Figure  7. 
Equating  the  horizontal  and  vertical  forces  to  zero  on  the  stress  block 
shown  yields  the  two  following  equations: 

A  A 

(o  )_  Ay  -  (t  )_  A*  +  F.  cos  0  As  -  F  sin  0  As  =  0 
M  xy  q  t  n 

A  A 

-  (0y)q  Ax  +  (xXy)g  Ay  +  Ft  sin  0  As  +  Fq  cos  0  As  =  0 

All  variables  refer  to  the  point  Q.  These  equations  can  easily  be 
siaqplified  to  the  following  form: 
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(16) 


2(tXy>Q  cos  0  sin  0  =  Fq 

(t  )a  (sin2  0  -  cos2  0)  =  F 

xy  Q  t 

Though  Equations  16  appear  to  contain  three  unknowns,  (^x)q>  (0y)q, 
and  (t  )q,  there  are  only  two  actual  unknowns,  since  by  substituting 
Equations  14  for  these  variables,  the  left-hand  sides  then  contain  only 
Pq  and  P  as  unknowns.  In  this  way  two  equations  in  two  unknowns  are 
written  for  each  point  Q  on  the  boundary  Tp. 

Likewise  at  each  point  Q  on  two  equations  in  the  unknowns  Pq 
and  P,.  are  obtained  since 

A 

U 

(17) 

V 

and  substitution  of  Equations  15  for  u^  and  Vq  leaves  the  left-hand 
sides  expressed  in  tens  of  PQ  and  Pfc. 

The  numerical  implementation  of  the  boundary  integral  equations 
proceeds  as  follows.  The  entire  boundary  is  subdivided  into  N  straight 
line  segments  at  the  center  of  which  the  pair  of  equations  in  either 
Equations  16  or  17  is  applied.  The  result  is  a  system  of  2N  boundary 
integral  equations  containing  the  unknown  tractions  PQ  and  Pt  acting 
over  each  segment.  P  and  Pfc  are  interpolated  at  the  same  N  points,  and 
2N  equations  in  2N  unknown  discrete  values  of  PQ  and  Pfc  are  obtained. 
Interpolation  of  the  unknown  traction  distribution  at  the  N  center 
points  is  consistent  with  assuming  the  tractions  to  be  constant  over  the 
segment.  Other  numerical  schemes  can  be  constructed  also.  For  example, 
the  unknown  tractions  can  be  interpolated  assuming  a  linear  or  parabolic 
distribution  over  the  segment.  In  this  study  the  constant  interpolation 
scheme  was  implemented,  although  the  equations  were  also  developed  for 
linear  interpolation. 

Once  the  interpolation  scheme  has  been  applied,  the  discrete  variables 

P  and  P,.  can  be  factored  outside  the  integrands.  The  remaining  integral 
n  t 

expressions,  when  evaluated  over  each  segment,  form  the  coefficients  of 


(ox)q  sin2  0  +  (oy)g  cos2  0  - 
[(oy)Q  -  (ox)q1  cos  6  sin  0  - 
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these  discrete  variables  in  the  system  of  algebraic  equations.  The 
integrations  were  carried  out  analytically  in  this  study,  as  was  the 
evaluation  of  the  singularity  contribution  given  in  Equations  14  and  15. 
The  integrals  could  have  been  evaluated  numerically  instead.  However, 
such  an  approach  essentially  would  have  resulted  in  a  numerical  integra¬ 
tion  over  each  segment  occurring  within  a  straight  line  approximation  of 
the  boundary  (i.e.,  an  approximation  within  an  approximation).  Again, 
the  details  of  the  analytical  integration  are  omitted.  The  system  of  2N 
linear  algebraic  equations  which  results  is 

a 

K  P  =  B  (18) 

SV  A<>  '  * 

The  square  coefficient  matrix  K  contains  the  results  of  the 
analytical  integration  of  the  singularities  in  2x2  blocks  on  the  diag¬ 
onal  and  the  analytical  integrations  of  the  terms  other  than  the  sin¬ 
gularities  in  the  off-diagonal  blocks.  The  vector  P  contains  unknown 

~  A 

variables  P  and  at  each  of  the  N  points.  The  vector  B  contains  the 
prescribed  boundary  values  of  either  the  tractions  (Fr  and  Ffc)  or  the 
displacements  (u  and  v)  at  each  of  the  N  points.  Some  details  of  the 
numerical  implementation  and  these  matrices  may  be  found  in  Appendix  A. 

The  solution  of  this  system  yields  2N  values  for  the  unknown  trac¬ 
tions  on  the  scribed  boundary  such  that  the  prescribed  boundary  condi¬ 
tions  are  imposed  in  the  infinite  sheet.  The  solution  is  written  as 

l  *  E'1  8  (is) 

In  this  study  the  solution  was  carried  out  by  a  Gaussian  elimination 
algorithm  with  partial  pivoting.  It  should  be  mentioned  that  for  homoge¬ 
neous  problems  the  coefficient  matrix  K  is  both  full  and  nonsymmetric. 

It  has  no  special  structure  that  could  otherwise  have  been  exploited  as 
in  finite  element  solutions  of  the  same  class  of  problems.  The  solution, 
Equation  19,  concludes  Step  2  of  the  indirect  boundary  element  procedure. 


18 


Step  3  is  described  as  follows.  The  prescribed  displacement  and 
stress  conditions  of  the  actual  problem  have  been  imposed  on  the  scribed 
boundary  in  the  infinite  sheet  of  the  auxiliary  problem.  The  Kirchoff 
uniqueness  theorem  is  invoked,  which  specifies  that  the  solution  to  the 
actual  problem  is  therefore  identical  to  the  solution  of  the  auxiliary 
problem.  As  a  result,  the  stresses  and  displacements  on  and  within  the 
scribed  boundary  of  the  infinite  sheet  are  identical  to  the  stresses  and 
displacements  on  and  within  the  finite  elastic  sheet  of  the  actual 
problem. 

Given  the  2N  values  for  Pq  and  P  ,  the  stress  or  displacement 
response  at  any  prescribed  field  point  Q  within  the  scribed  boundary  can 
be  found  using  the  appropriate  Green's  functions,  Equations  10  or  11. 
Referring  to  Figure  8,  the  stress  and  displacement  responses  at  a  pre¬ 
scribed  field  point  Q  are  given  by  the  following  equations.  The  stresses 
at  Q  are 
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The  displaceaents  at  Q  are 
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CWIPUTATIONAL  PERFORMANCE  OF  THE  INDIRECT  BOUNDARY  ELEMENT  METHOD 


Beaai  Resting  on  an  Elastic  Foundation 

A  numerical  study  was  conducted  by  comparing  computed  results  from 
two  indirect  BEM  programs  with  theoretical  solutions  for  beam-on-an- 
elastic-foundation  problems.  The  first  program  was  written  in  FORTRAN 
and  implemented  on  a  Cyber  175  mainframe  computer.  The  second  program, 
developed  independent  of  the  first,  was  written  in  Applesoft  BASIC  on  an 
Apple  II  Plus  microcomputer.  Other  boundary  element  computer  programs 
have  been  written  for  microcomputers  (Ref  14  and  15),  but  the  authors 
are  unaware  of  any  based  upon  the  indirect  method  for  stress  analysis. 

There  are  soaie  differences  in  these  two  programs.  The  FORTRAN  ver¬ 
sion  contains  closed-form  solutions  for  full-span  uniform  and  bilinear 
varying  loading,  requires  batch  input,  and  assumes  free-end  boundary 
conditions.  The  BASIC  version  approximates  both  uniform  and  linear 
varying  loads  (partial  or  full  span)  with  equally  spaced  concentrated 
loads.  The  input  is  interactive,  and  each  set  of  boundary  conditions 
can  be  specified  in  terms  of  displacement,  rotation,  moment,  or  shear. 
The  larger  real  word  size  of  the  Cyber  175  computer  (60  bits  versus 
40<bit8  for  the  Apple)  and  the  closed-form  solution  for  linear  varying 


loads  contribute  to  the  better  accuracy  of  FORTRAN  versions.  The  BASIC 
version  has  more  flexibility  with  interactive  input,  more  general  loading, 
and  the  ability  to  model  more  complex  boundary  conditions. 

The  problems  shown  in  Figures  9  through  14  are  the  subjects  of  the 
numerical  study  for  one-dimensional  applications.  These  problems  include 
concentrated  loads,  a  partial  uniform  load,  and  exploit  symmetry  condi¬ 
tions.  Appendix  B  includes  additional  problems  for  the  BASIC  version, 
showing  linear  varying  loads  and  asymmetry  modeling.  The  listing  of  the 
BASIC  program  is  included  in  Appendix  C. 

As  shown  in  Figure  9,  the  first  problem  consist  of  a  50-kip  load  on 
the  center  of  a  beam.  The  moment  and  deflection  at  the  center  of  the 
beam  are  compared  with  theory.  The  results  show  both  programs  agreeing 
exactly  (to  the  indicated  accuracy)  with  theory  for  deflection  at  the 
midspan.  Both  programs  give  accurate  results  for  the  moment  value,  the 
FORTRAN  version  being  the  most  accurate.  Figure  10  demonstrates  the  use 
of  symmetry  in  the  solution  of  the  same  problem. 

The  mainframe  version  of  the  program  executed  100  times  faster  than 
the  microcomputer  version.  Using  a  compiled  version  of  the  BASIC  program 
nearly  doubled  its  speed.  Though  the  programs  were  constructed  indepen¬ 
dently,  this  does  give  a  rough  estimate  of  the  speed  difference  between 
the  mainframe  and  this  generation  of  microcomputer. 

Figure  11  illustrates  the  second  example  involving  two  symmetric 
concentrated  loads.  The  moment  and  deflection  at  the  center  of  the  beam 
are  compared  with  theory.  Figure  12  shows  the  model  of  the  problem 
incorporating  symmetry  boundary  conditions.  Both  programs  give  equally 
good  results.  The  symmetric  modeling  gives  the  same  values  as  the  full 
beam  model. 

Figures  13  and  14  illustrate  the  third  example  involving  a  partial 
uniform  load  symmetric  about  the  centerline.  The  moment  and  deflection 
at  the  center  of  the  beam  and  the  slope  at  the  ends  of  the  beam  are 
compared  with  theory.  The  BASIC  program  gives  better  results  for  the 
end  slopes,  while  the  FORTRAN  program  gives  better  results  for  the 
centerline  moment.  The  accuracy  of  the  bending  moment  calculation  from 
the  BASIC  version  using  symmetry  is  equal  to  that  from  the  FORTRAN 


version.  The  full  model  uses  20  concentrated  loads  to  approximate  the 
partial  uniform  load.  The  symmetry  model  uses  both  10  and  20  concen¬ 
trated  loads  to  approximate  the  partial  uniform  load.  All  models  indi¬ 
cate  about  the  same  accuracy,  indicating  that  10  concentrated  loads  are 
adequate  to  represent  the  partial  uniform  load  in  the  symmetry  model. 

The  interactive  mode,  allowable  with  the  BASIC  version,  makes  the  micro¬ 
computer  version  more  valuable  as  an  engineering  tool. 

Other  examples  presented  in  Appendix  C  include  a  partial  span, 
linear  varying  load,  and  antisymmetric  loading.  The  inversion  test  (see 
Appendix  B  for  an  explanation)  indicated  that  the  solution  might  have 
been  inaccurate.  A  possible  problem  with  the  boundary  integral  method 
is  that  the  principal  coefficient  matrix  ([G]  in  this  case)  contains 
terms  that  can  vary  by  orders  of  magnitude.  The  printout  of  the  [G] 
matrix  indicates  zero  terms  (to  three  decimal  places)  for  the  displacement 
and  rotation  rows.  The  matrix  might  require  preconditioning  to  ensure 
consistently  good  numerical  accuracy. 

Of  the  11  problems  studied,  all  but  one  indicate  the  BEM  to  be  very 
accurate  for  the  beam-on-an-elastic-foundation  problems.  The  only 
approximations  made  in  the  algorithms  were  the  reduction  of  continuous 
loads  to  a  series  of  concentrated  loads  in  the  BASIC  program. 

Two-Dimensional  Elastostatics 

A  brief  numerical  study  was  conducted  on  three  simple  example  prob¬ 
lems  in  which  the  solutions  from  an  indirect  BEM  program,  a  direct  BEM 
program,  and  theory  were  compared.  Both  programs  incorporated  constant 
boundary  elements.  The  indirect  BEM  program  (BIM2D)  was  developed  as 
part  of  this  study,  while  the  direct  BEM  program  (PGM18)  was  developed 
by  C.  A.  Brebbia  at  Southampton  and  is  described  in  Reference  16. 

Appendix  D  contains  a  listing  of  the  BIM2D  program. 

In  the  boundary  element  method,  most  of  the  computer  time  is  spent 
integrating  the  boundary  integrals  over  the  elements  to  form  the  system 
matrix  K  and  calculating  internal  responses  depending  on  how  much  inter- 

A/ 

nal  information  is  sought.  In  the  finite  element  method  the  solution  of 
the  system  of  algebraic  equations  is  usually  the  longest  calculation. 
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In  addition  to  accuracy  considerations,  BIM2D  was  developed  using  ana¬ 
lytical  integrations  over  the  elements  to  minimize  the  cost  of  boundary 
integrations.  This  offered  better  accuracy  for  a  given  cost.  Integra¬ 
tions  within  the  direct  BEM  program  are  carried  out  numerically  except 
for  the  singularity  coefficients  -  diagonal  blocks  of  the  system  matrix. 
Several  observations  drawn  from  the  numerical  study  are  discussed  below. 

The  ripple  effect  (as  shown  in  Figures  15  and  16)  occurs  near  a 
boundary  where  a  relatively  few  number  of  elements  are  used.  The  greatest 
deviations  occurred  near  abrupt  geometry  discontinuities  (corners)  and 
loading  discontinuities . 

The  first  problem  studied  was  a  1-inch-thick  square  plate  in  hydro¬ 
static  tension.  As  shown  in  Figure  17  the  model  consisted  of  12  ele¬ 
ments  of  unit  length  on  a  side,  all  with  a  prescribed  normal  traction  of 
1,000  psi.  Internal  stresses  were  calculated  in  the  upper  quarter  of 
the  plate.  The  stresses  at  four  lines  of  30  response  points  each  are 
shown  plotted  in  Figure  15. 

Except  for  the  edge  region,  within  one  element  length  of  the  bound¬ 
ary,  the  direct  method  algorithm  is  more  accurate.  Excluding  the  edge 
region,  the  indirect  algorithm's  error,  5.6%  and  less,  is  on  the  order 
of  10  times  that  of  the  direct  algorithm's  error.  In  the  edge  region 
the  indirect  algorithm  is  more  accurate,  and  this  is  the  critical  region 
in  most  applications.  The  direct  method  behaves  pathologically  near  the 
boundary,  while  the  indirect  algorithm  is  relatively  stable.  Both 
methods  exhibit  a  ripple  effect  or  artificial  stress  oscillation,  but  it 
is  much  more  pronounced  in  the  direct  algorithm.  It  is  interesting  that 
the  normal  stress  computed  from  the  direct  method  changes  direction  of 
instability  near  the  corner.  That  is,  a  goes  one  way  and  a  goes  the 

y  * 

other . 

The  indirect  method  always  underestimates  the  response  where  the 
stresses  are  largest.  For  most  response  lines  examined  the  response 
begins  with  a  nearly  constant  accuracy,  decreases  in  accuracy,  and  then 
recovers  near  the  boundary.  The  presence  of  the  boundary  influences  a 
larger  region  in  the  indirect  algorithm,  but  this  effect  is  less  severe. 
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The  second  problem  studied  was  a  1-inch-thick  rectangular  plate 
subjected  to  1,000  psi  hydrostatic  tension.  As  shown  in  Figure  18,  the 
model  consisted  of  two  20-element  sides  and  two  4-element  sides.  Each 
element  had  a  length  of  unity.  The  internal  responses  were  calculated 
in  the  upper  quarter  of  the  plate.  The  stresses  along  four  lines  of 
50  response  points  each  were  calculated  and  are  shown  plotted  in  Figure  16 
This  problem  was  also  the  subject  of  a  numerical  study  conducted  by  John 
Bode  (Ref  12). 

In  this  example  the  computed  stress  response  was  sampled  within 
one-tenth  of  an  element  length  of  the  plate  edge.  The  resulting  ripple 
effect  is  much  more  pronounced  than  it  was  in  the  square  plate  results. 

The  indirect  method  exhibits  a  ripple  effect  of  about  ±1.0%.  This  is 
overshadowed  by  the  strong  divergence  from  the  solution  exhibited  by  the 
direct  method  algorithm.  As  with  the  first  example  the  direct  method  is 
more  accurate  excluding  the  edge  region.  However,  a  high  boundary/domain 
ratio  problem  (in  two  and  three  dimensions)  requires  an  algorithm  that 
is  more  accurate  in  the  edge  region,  since  this  region  comprises  a 
greater  portion  of  the  domain.  The  indirect  method  continues  to  exhibit 
a  loss  and  then  recovery  in  accuracy  near  the  boundary  except  in  the 
corner  region. 

The  corner  region  is  the  area  of  strongest  divergence  for  both 
methods.  Since  the  term  "constant  element"  alludes  to  the  artificial 
boundary  loads,  the  BEM  will  not  give  an  exact  solution  for  a  uniform 
stress  field  (unlike  the  finite  element  method  with  a  constant  strain 
element,  for  example).  In  the  first  two  problems  considered  the  pre¬ 
scribed  boundary  conditions  were  constant  (hydrostatic  tension) .  The 
indirect  BEM  solves  the  problem  in  the  infinite  domain,  and  the  gradient 
of  the  artificial  boundary  stresses  is  noted  to  be  high  near  the  cor¬ 
ners.  It  is  believed  that  the  artificial  boundary  stresses  tended  to 
increase  near  the  corners  to  maintain  the  square  shape  of  the  plate. 

A  concentration  of  elements  near  the  corners  yields  a  more  accurate 
solution  by  better  modeling  of  the  artificial  tractions.  However,  for  a 
model  with  a  fixed  number  of  elements,  overconcentration  near  the  corners 
increases  the  ripple  effect  near  the  large  elements.  One  method  for 


refining  a  boundary  mesh  might  be  based  on  maintaining  a  constant  value 
for  the  total  artificial  traction  per  element.  Bode's  study  showed  that 
beveling  the  corners  was  not  effective  with  constant  elements. 

In  Bode's  study,  96  well-placed  elements  gave  average  values  almost 
as  good  as  288  equally  spaced  elements.  The  ripple  effect,  however,  was 
more  prominent.  Central  processing  unit  (CPU)  time  is  the  major  factor 
in  the  number  of  elements.  The  following  table  shows  the  relationship 
between  the  number  of  elements  and  CPUs  in  Bode's  study.  The  third 
problem  in  our  study,  a  stress  concentration  problem,  did  not  show  cost 
as  strongly  dependent  on  the  number  of  elements  (see  Figure  19). 


Number  of  Elements 

Normalized  CPUs 

96 

1.0 

192 

3.5 

288 

8.65 

Forty-eight  elements  of  unit  length  were  used  in  both  the  square 
and  rectangular  plate  problems.  At  the  response  lines  half  an  element 
length  away  from  the  edge,  the  results  from  the  square  plate  problem 
were  more  accurate.  At  the  center  of  the  plate,  stress  errors  of  2%  to 
3%  were  observed  for  the  square  plate  and  5%  to  7%  were  observed  for  the 
rectangular  plate.  Thus,  the  proximity  of  edges  has  an  adverse  effect 
on  stress  response  accuracy. 

The  higher  inaccuracies  in  the  rectangular  plate  problem  are  prob¬ 
ably  due  to  each  response  point  being  relatively  closer  to  more  edges 
and  a  greater  number  of  elements.  In  this  example  the  BEM  slightly 
decreased  in  accuracy  as  the  narrowness  of  the  domain  increased.  For  a 
given  accuracy  the  difference  in  analysis  cost  for  the  two  plate  problems 
would  not  vary  much  with  the  BEN.  However,  the  cost  would  be  expected 
to  vary  if  the  finite  element  method  (FEM)  were  to  be  used. 

The  third  problem  studied  was  a  1-inch-thick  square  plate  with  a 
circular  cut  out  that  involves  high  stress  gradients  or  stress  concen¬ 
trations.  Uniform  compression  stress  of  1  psi  was  applied  to  two  opposite 
edges.  Various  numbers  of  elements  were  used  to  model  the  external  and 
internal  boundaries,  but  in  all  cases  the  elements  for  a  given  boundary 


were  of  uniform  size.  The  various  models  are  indicated  in  Figure  19.  A 
vertical  and  a  horizontal  line,  of  37  response  points  each,  extended 
from  the  hole  to  the  outer  boundary.  Figure  20  shows  the  stress  response 
for  the  different  models.  The  theoretical  solution  shown  assumes  the 
elastic  plate  to  be  of  infinite  extent. 

Since  the  response  lines  did  not  closely  parallel  an  edge,  no 
ripple  effect  was  observed.  As  with  the  other  two  problems  the  direct 
algorithm  diverges  strongly  at  about  half  an  element  length  from  the 
boundary.  Because  of  the  strong  divergence,  the  direct  algorithm  was 
only  used  for  comparison  with  the  coarse  model,  model  1.  For  the  indirect 
algorithm  the  responses  are  relatively  accurate,  diverging  slightly  in 
the  edge  regions.  The  models  with  the  finer  element  subdivision  on  a 
given  boundary  tend  to  give  the  most  accurate  results  near  that  boundary. 
The  level  of  subdivision  on  distant  boundaries  has  only  a  secondary 
effect.  Figure  20b  illustrates  the  variation  of  ax  along  the  vertical 
response  line.  For  this  response  the  localized  effect  of  the  boundaries 
is  evident.  Near  the  circular  edge  where  the  stress  gradient  is  highest, 
models  3  and  4  show  the  best  results  and  contain  approximately  4%  error. 
Models  1  and  2  have  an  error  of  between  5%  and  6%.  The  number  of  elements 
used  on  the  square  has  only  a  secondary  effect. 

Models  2  and  4  show  the  best  results  along  the  square  edge  and 
contain  approximately  8%  error.  Models  1  and  3  have  an  error  of  approx¬ 
imately  13%.  Again,  the  number  of  elements  used  to  model  the  distant 
boundary,  in  this  case  the  circular  edge,  has  only  a  secondary  effect. 

Figure  20c  illustrates  the  variation  of  along  the  vertical 
response  line.  In  theory  0^  should  be  zero  at  both  boundaries.  As 
awntioned  earlier  the  theoretical  curve  corresponds  to  an  infinite 
plate,  and  this  is  the  main  source  of  discrepancy.  For  models  1  and  3 
(both  having  a  28-element  square  boundary)  a  slight  decrease  in  accuracy 
with  a  recovery  near  the  square  boundary  is  exhibited  again. 

Figure  20d  illustrates  the  variation  of  0x  along  the  horizontal 
response  line.  In  theory  this  stress  should  be  0  at  the  circle  and 
-1  psi  (the  applied  traction)  at  the  outer  boundary.  The  theoretical 
solution  is  approaching  -1  psi  slower  than  the  BEM  solution. 


Figure  20e  illustrates  the  variation  of  along  the  horizontal 
response  line.  The  accuracy  trends  previously  mentioned  are  slightly 
deviated  from  near  the  circle  because  in  this  example  model  2  gives 
slightly  better  results  than  model  3  while  having  fewer  elements  on  the 
circular  boundary.  The  maximum  error  for  a on  the  circular  boundary 
was  20%  for  this  response  line.  But  it  must  be  noted  that  is  of  les¬ 
ser  importance  in  the  stress  concentration  example. 

In  general,  the  constant-element  BEM  exhibited  a  good  ability  to 
model  stress  gradients.  Unlike  the  uniform  stress  field  problems,  the 
stress  gradient  problem  showed  the  BEM  to  overestimate  stress  values 
(i.e.,  to  be  conservative).  The  boundary  effects  are  observed  to  be 
local  and  relative  to  the  element  size.  Distant  boundaries  are  seen  to 
have  only  a  secondary  effect. 

SUMMARY  AND  CONCLUSIONS 

The  indirect  boundary  element  method  has  been  investigated  for 
application  to  one-  and  two-dimensional  elastostatics  problems  in  struc¬ 
tural  analysis.  The  theoretical  basis  of  the  method  has  been  described 
by  beginning  with  the  simple  problem  of  a  beam  resting  on  an  elastic 
foundation  and  then  extending  the  theory  to  a  more  useful  class  of 
problems  pertaining  to  the  plane  stress  and  plane  strain  analysis  in 
elastostatics . 

The  numerical  implementation  of  the  theoretical  formulation  was 
illustrated  with  the  development  of  computer  programs  for  both  small  and 
large  computers.  Stress  analysis  capability  using  these  computer  pro¬ 
grams  was  assessed  by  comparison  of  results  to  theoretical  solutions  and 
to  results  from  another  computer  program  that  is  based  on  the  alterna¬ 
tive  direct  boundary  element  formulation. 

Three  basic  concepts  constitute  the  theoretical  formulation  of  the 
indirect  boundary  element  method:  Green’s  functions,  superposition,  and 
the  Kirchhoff  uniqueness  theorem. 


Green's  functions  are  akin  to  the  more  familiar  influence  functions, 
and  they  give  the  stress  and  displacement  response  at  an  arbitrary  field 
point  due  to  a  unit  force  at  an  arbitrary  source  point  in  the  domain. 
These  functions  are  classical  results  from  the  theory  of  elasticity,  and 
those  which  were  used  herein  apply  to  infinite  domains  only.  They  sug¬ 
gest  the  reformulation  of  actual  problems  in  the  finite  domain  in  terms 
of  problems  in  the  infinite  domain.  This  gives  rise  to  the  auxiliary 
problem,  which  is  solved  in  place  of  the  actual  problem,  and  thus  the 
origin  of  the  word  "indirect"  in  the  indirect  boundary  element  method. 

Through  superposition  of  the  Green's  functions,  a  set  of  boundary 
integral  equations  is  constructed  whose  solution  imposes  the  prescribed 
boundary  conditions  on  a  contour  in  the  infinite  domain  that  is  identi¬ 
cal  to  the  boundary  contour  of  the  actual  finite  problem. 

Kirchhoff's  uniqueness  theorem  then  requires  that  the  solution  of 
the  auxiliary  problem  be  identical  to  the  solution  of  the  actual  problem. 
The  solution  of  the  auxiliary  problem  is  accomplished  numerically.  The 
boundary  is  discretized  into  N  straight  line  segments  over  which  unknown 
artificial  normal  and  tangential  stress  tractions  are  interpolated.  In 
this  study  a  constant  value  for  the  unknowns  was  assumed  over  each 
segment.  The  integration  of  the  boundary  integrals  over  a  segment  was 
carried  out  analytically,  and  the  results  are  the  coefficients  of  the  2N 
unknown  tractions.  A  linear  system  of  2N  algebraic  equations  in  the  2N 
unknowns  occurs  and  is  solved  by  Gaussian  elimination.  The  coefficient 
matrix  has  little  exploitable  structure  for  the  homogeneous  examples 
considered.  It  is  both  fully  populated  and  nonsymmetric.  -Once  the 
artificial  tractions  along  the  boundary  are  known,  the  stress  and  dis¬ 
placement  response  may  be  determined  by  superposition. 

The  results  from  the  rather  simple  two-dimensional  numerical  exper¬ 
iments  carried  out  in  this  study  suggest  that  the  boundary  element 
methods  are  susceptible  to  accuracy  deterioration  within  one  element 
length  of  the  boundary  or  edge  of  the  domain.  However,  within  the 
domain  the  accuracy  is  very  satisfactory,  usually  always  within  5%  of 
exact  solutions. 


Very  near  the  edges,  however,  a  ripple  effect  occurs  in  the  computed 
stress  values.  The  ripple  effect  is  aa  oscillation  in  the  stress  magni¬ 
tude  about  the  exact  value.  The  oscillation  along  a  line  paralleling  an 
edge  increases  as  a  corner  of  the  domain  is  approached.  That  is,  as 
another  edge  is  approached,  the  ripple  effect  becomes  worse. 

Since  the  critical  region  of  the  ripple  effect  is  within  one  element 
length  of  an  edge,  element  size  is  a  factor  in  achieving  accuracy. 

Local  element  size  gradation  can  be  employed  to  achieve  improved  accuracy 
near  an  edge,  while  distant  element  size  has  little  influence  on  the 
accuracy  in  the  vicinity  of  this  edge. 

Inasmuch  as  significant  stresses  often  occur  on  boundaries,  fine 
subdivisions  of  those  boundaries  may  be  necessary  to  achieve  desired 
accuracy. 

With  regard  to  edge  effects,  the  indirect  boundary  element  method 
faired  better  than  the  particular  direct  boundary  element  implementation 
that  was  available  for  comparison. 

The  results  pertaining  to  the  accuracy  of  the  indirect  boundary 
element  method  to  correctly  capture  stress  gradients  were  very  encourag¬ 
ing,  even  for  the  constant  stress  elements  employed.  The  results  suggest 
that  the  method  may  be  a  very  economical  analysis  tool  for  determining 
stress  concentration  factors  in  elastostatics. 

Compared  to  finite  element  swthods,  the  necessary  input  data  require¬ 
ments  are  less,  and  smaller  matrices  result.  The  boundary  element 
methods  are  therefore  more  suited  to  small  computers.  This  would  then 
allow  the  personal  interactive  advantages  of  microcomputers  to  further 
enhance  the  boundary  element  methods  as  effective  stress  analysis  tools. 

RECOMMENDATION 

From  the  results  of  this  study,  it  is  believed  that  a  combined 
finite  element  and  boundary  element  computer  program  may  prove  successful 
in  contributing  to  the  reduction  in  the  high  costs  now  associated  with 
nonlinear  finite  element  programs.  A  two-dimensional  program  could  be 


developed  and  used  to  evaluate  potential  advantages  in  the  solution  of 
nonlinear  problems,  particularly  soil-structure  interaction  problems 
where,  for  example,  a  buried  structure  along  with  some  surrounding  soil 
may  be  behaving  nonlinearly  and  the  remaining  half-space  soil  is  behaving 
linearly.  However,  the  proper  theoretical  and  numerical  treatment  of 
the  interface  equations  must  be  thoroughly  evaluated  and  then  implemented. 
This  step  is  a  necessary  prerequisite  to  an  effective  implementation. 

The  recommendation  is  to  undertake  the  theoretical  and  numerical  treatment 
of  coupling  the  finite  element  and  indirect  boundary  element  methods. 
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Figure  6.  Actual  two-dimensional  region  embedded  in 
infinite  two-dimensional  region  -  the 
auxiliary  problem. 
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Figure  9.  Full  model  of  a  concentrated  load  on  centerline. 
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Figure  10.  Symmetric  model  of  a  concentrated  load  on  centerline. 
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Figure  11.  Full  model  of  symmetric  concentrated  load. 
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Figure  12.  Symmetric  model  of  symmetric  concentrated  loads. 
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Figure  13.  Full  model  of  partial  uniform  load. 
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Figure  14.  Symmetric  model  of  partial  uniform  load. 
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Figure  15b.  Normal  stress,  oy,  at  tour  response  lines. 


Figure  16a.  Rectangular  plate  model,  normal  stress,  a. 
at  four  response  lines. 


Figure  16c.  Shear  stress,  f—,,  at  tour  response  lines. 


Figure  18.  Rectangular  plate  model. 


0  Indirect 


Figure  20d.  ffx  along  horizontal  response  line. 


Figure  20e.  oy  along  horizontal  response  line. 


Appendix  A 


ASPECTS  OF  THE  NUMERICAL  IMPLEMENTATION 

This  appendix  includes  further  detail  on  the  numerical  implementa¬ 
tion  of  the  indirect  boundary  element  method  (BEM).  Each  section  explains 
one  aspect  that  was  omitted  for  clarity  in  the  text  of  the  report  on 
two-dimensional  elastostatics.  Though  some  of  the  detail  is  unique  to 
the  BIM2D  program,  the  overall  methodology  is  general  in  nature. 

SIMPLIFIED  BOUNDARY  INTEGRAL  EQUATIONS 

Equation  14a  is  rewritten  as  follows.  The  contribution  of  the 
singular  element  is  written  as  a  product  of  influence  coefficients  and 
discrete  artificial  tractions.  The  integration  along  the  remainder  of 
the  boundary  is  written  as  the  summation  of  the  integrations  along  the 
nonsingular  elements.  In  the  present  study,  these  integrations  were 
carried  out  analytically,  but  these  details  are  omitted  and  only  the 
results  are  given.  Using  the  following  substitutions. 


A-l 


Equation  14a  can  be  rewritten  in  the  following  form: 


<°,>Q  ■  SI“Q,Q  ♦  SXIQ,Q  % 


N 

-  E 

p=i 

P*Q 


(-  P  sin  0  +  P  cos  6)  i 


+  (Pnp  cos  e  yt?  sin  6)  Y2Qip 
Combining  terms  with  respect  to  the  boundary  tractions  gives 

<VQ  ■  SalQ,Q  \  *  SXTQ,Q  PtQ 


N 

-  E 

p=i 

P*Q 


Pnp(-sin  0  YlQp  +  cos  0  y2q>p 


+  P^(cos  y1q>p  +  sin  8  Y2Q>p) 


(A-3) 


(A-4) 


The  integration  along  the  remainder  of  boundary  can  now  be  written  in 
coefficient  form.  The  nonsingular  influence  coefficients  are  given  by: 


sxnq,p  =  '  8in  0  y1q,p  *  cos  y2q,p 
sxtq,p  =  cos  y1q,p  +  8in  y2q,p 

Thus,  Equation  A-4  is  further  simplified  to: 

N  N 

(°Jn  =  E  K  SXNo  p  +  E  pt  SXTo  p 
*  Q  ft  “p  Q*P  ft  lP  Q»p 


(A-5a) 

(A-5b) 


(A-6a) 


where  the  singular  coefficients,  subscripted  Q,Q,  are  given  by  Equa¬ 
tions  A-l,  and  the  nonsingular  coefficients,  subscripted  Q,P,  are  given 
by  Equations  A-S. 


Io  a  siailar  Banner  (ay)q»  ^xy'Q*  (u)q»  and  (v'q  can  be  rewritten 


<V«  =  £  %  smq,p  *  I,  %  S¥Tq,p 


(A-6b) 


<  Vo  *  E>.  "'I',  11 

*y  ^  p=i  np  p=i  cp 


(A-6c) 


(u)n  =  Z  P  UN  +  £  Pf  ™ 
Q  W  “P  £l  fcp 


(A-6d) 


WQ  *  E  K  w  ♦  E  rt  ” 

W  P=1  P  P=1  P 


(A-6e) 


The  coefficients  SYN_  p  and  SYT_  _  are  defined  in  a  manner  similar  to 
that  shown  in  Equations  A-l  and  A-5  using  Equations  14b;  the  coeffi¬ 
cients  TN  and  TT  use  Equation  14,  the  coefficients  UN  and  UT  use  Equa¬ 
tion  15a,  and  the  coefficients  VN  and  VT  use  Equation  15b. 


PRESCRIBED  BOUNDARY  TRACTIONS  IN  TERMS  OF  THE  ARTIFICIAL  BOUNDARY 
TRACTIONS 

Equation  16  gives  the  prescribed  boundary  tractions  in  terms  of  the 
unknown  stresses  (°x)q>  ^ay^Q>  an<*  ^xy^Q'  Writing  the  unknown  stresses 
in  terms  of  the  artificial  boundary  tractions,  using  Equations  A-6a,  b, 
and  c,  gives: 


£  P  [sin2  B(SXN.  „)  +  cos2-  8(SYWn  -)  -  2  sin  6  cos  8  (TH.  _)) 

“p  Q*P  Q>P  Q’F 

♦  P^Isin2  0(SXTq  p)  +  cos2  8(SYTq  -  2  sin  8  cos  8  (TNQ>p)]  = 


(A-7a) 


Jp^J-cos  8  sin  8  (SYNq  p)  +  cos  0  sia  8  (SYNq  p)  +  (cos2  0  -  sin2  0)(THq  p)) 
♦  P,  [-cos  e  sin  6  (SXT„  «)  ♦  cos  8  sin  0  (SYT-  „)  +  (cos2  0  -  sin2  8)(TT-  p))l 


=  Ft  (A-7b) 


The  quantities  within  the  square  brackets  are  the  coefficients  of  the 
matrix  K.  Thus,  Equations  A- 7  can  be  rewritten  as: 


(A-8a) 

(A-8b) 


In  this  forsi  it  is  apparent  that  the  influence  of  the  artificial  boundary 
tractions  at  each  element  P  with  the  prescribed  boundary  tractions  at 
element  Q  results  in  a  block  of  four  teras  (a  2x2  subaatrix)  in  the 
coefficient  aatrix. 


PRESCRIBED  BOUNDARY  DISPLACEMENTS  IN  TERMS  OF  THE  ARTIFICIAL  BOUNDARY 
TRACTIONS 


Equation  17  gives  the  prescribed  boundary  displacements  (u,  v)  in 
teras  of  the  unknown  boundary  displaceaents  (uq,  Vq).  Writing  these 
displaceaents  in  teras  of  the  artificial  boundary  tractions,  using 
Equations  A-6d  and  e,  gives: 


(A-9a) 

(A-9b) 


The  quantities  within  the  square  brackets  are  the  coefficients  of  the 
aatrix  K.  Thus,  Equations  A-9  can  be  rewritten  in  the  same  form  as 
Equations  A-8: 


A-4 


I 


N 

where  i  =  2xQ  -  1 
j  =  2xP  -  1 


tp  x+l,j+lj 


Again,  the  interaction  of  an  ordered  pair  of  elements  results  in  a  block 
of  four  terms  in  the  coefficient  matrix. 


SYSTEM  OF  EQUATIONS 

Equations  A-8  and  Equations  A- 10  each  yield  a  set  of  two  equations 
in  2N  unknowns.  Either  set  of  equations  can  be  written  in  matrix  notation 


K.  .  K.  .  K.  ,  K.  . 

1,1  1,2  i,3  1,4 

Ki+l,l  Ki+1,2  Ki+1,3  Ki+1,4  * 


*  Ki,2N-l  Ki,2N  I  “l 


•  Ki+1 ,2N-1  Ki+1,2N  II  ll 


(A- 11) 


where  i 


=  2xQ  -  1 


=  prescribed  x  displacement  or  normal  traction  at 
element  Q 


=  prescribed  y  displacement  or  tangential  traction 
at  element  Q 


For  each  straight-line  element  around  the  model  boundary,  either 
the  boundary  displacements  or  tractions  are  prescribed.  Thus,  each 
«l«*c&t  produces  two  equations.  The  total  system  of  equations  (Equation  18) 
is  formed  by  combining  N  pairs  of  equations,  such  as  Equations  A-10. 

Thus,  2N  unknown  artificial  boundary  tractions  are  expressed  in  terms  of 
2H  prescribed  boundary  tractions  or  displacements. 
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BIM1D  DOCUMENTATION  AND  LISTING, 
BASIC  VERSION 


PURPOSE 

This  appendix  illustrates  how  the  nature  of  the  boundary  integral 
■ethod  lends  itself  to  the  solution  of  structural  analysis  problems  on  a 
Microcomputer. 


SCOPE 


The  BIM1D  program  solves  the  problem  of  a  beam  on  an  elastic  founda¬ 
tion.  The  loading  is  limited  to  concentrated  loads;  however,  it  is  set 
up  to  generate  loads  to  approximate  linear  varying  continuous  loads. 

The  boundary  condition  at  each  end  of  the  finite  beam  can  be  specified 
by  entering  the  value  of  two  of  the  following  quantities: 

•  shear 

•  moment 

•  displacement 

•  rotation 

This  allows  modeling  half  of  the  beam  for  symmetrical  and  asymetrical 
loading  cases  (see  Appendix  C). 


METHODOLOGY 

The  theory  for  the  program  is  explained  in  detail  in  the  body  of 
the  report.  The  following  discussion  briefly  explains  the  major  sub¬ 
routines  . 

Input  Routine  (Beginning  Statement  Number  1050) 

The  input  routine  is  interactive  only  to  the  extent  that  it  prompts 
the  user  for  input.  It  does  not  check  the  range  of  the  user's  input 
(with  the  exception  of  allowing  a  maximum  of  20  uses  of  defined  concen¬ 
trated  loads)  and  does  not  allow  editing  of  input. 

Generator  Routine  (Beginning  Statement  Number  920) 

The  generator  routine  generates  concentrated  loads  to  approximate 
the  continuous  loads.  These  loads  are  added  to  the  concentrated  load 
arrays  (both  positions  and  values). 
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Matrices  Development  Routine  (Beginning  Statenent  Number  600) 


The  matrices  development  routine  initially  adjusts  the  position  of 
all  concentrated  loads  such  that  no  load  is  within  a  given  distance, 
epsilon,  from  either  boundary.  The  routine  then  develops  both  the  [G] 
and  (Hp)  matrices. 

j?  Matrix  Inversion  Routine  (Beginning  Statement  Number  750) 

The  matrix  inversion  routine  inverts  the  4x4  [G]  matrix  and  per¬ 
forms  an  accuracy  test.  In  the  accuracy  test  the  flexibility  matrix  [G] 
is  multiplied  by  its  inverse  to  give  the  identity  matrix.  The  sum  of 
all  terms  in  the  identity  matrix  is  printed  on  the  screen  and  should 
approach  4  (the  order  of  the  matrix).  If  [G]  is  a  singular  matrix  the 
program  ends. 

Solve  Routine  (Beginning  Statement  Number  730) 

The  solve  routine  solves  for  the  boundary  forces  and  moments  such 
that  when  they  are  applied  to  the  infinite  beam,  the  points  within  the 
boundaries  respond  as  if  they  were  on  the  finite  beam. 

Response  Input  Routine  (Beginning  Statement  Number  1800) 

The  response  input  routine,  as  the  initial  input  routine,  is  inter¬ 
active  only  to  the  extent  of  prompting  the  user  for  input.  It  allows  a 
maximum  of  50  response  points  and  does  not  allow  a  response  point  outside 
of  the  span. 

Response  Routine  (Beginning  Statement  Number  250) 

The  response  routine  initially  adds  the  boundary  moments  to  the 
applied  moment  array  (it  is  set  up  in  this  manner  so  that  applied  mosients 
are  a  simpler  addition);  similarly,  the  boundary  forces  are  added  to  the 
applied  force  array.  Each  response  point  is  then  compared  to  each 
concentrated  load  position.  If  they  coincide  the  response  point  is 
moved  to  the  right  by  the  amount  epsilon  (response  points  at  L  are  moved 
to  the  left).  With  singularities  avoided,  the  responses  can  be  cal¬ 
culated.  For  each  response  point,  (kB)  and  [kM]  are  developed  then 
multiplied  by  (B)  and  (M),  respectively.  This  routine  prints  on  the 
screen  the  response  point  that  is  being  calculated. 

Ouput  Routine  (Beginning  Stateatent  Number  3000) 

Tile  output  routine  formats  and  prints  the  calculated  results  and 
user  input.  A  "DEBUG  PRINTOUT"  is  optional  as  explained  in  the  User 
Instructions.  Very  large  numbers  can  cause  an  illegal  quantity  error. 
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USER  INSTRUCTIONS 


I.  PROGRAM  INPUT 

A.  Model  Input. 

1.  Enter  the  foundation  spring  constant. 

2.  Enter  the  modulus  of  elasticity  of  the  beam. 

3.  Enter  the  moment  of  inertia  of  the  beam. 

4.  Enter  the  span  of  the  beam. 

5.  Enter  the  number  of  continuous  loads. 

For  Each  Continuous  Load 


5a.  Enter  the  position  and  value  of  the  left  end  of  the 
continuous  load  in  the  format  position,  value. 

5b.  Repeat  Step  6  for  the  right  end  of  the  continuous  load. 
5c.  Enter  the  number  of  concentrated  loads  to  approximate  the 
continuous  load. 

6.  Enter  the  number  of  concentrated  loads. 

For  Each  Concentrated  Load 


Enter  the  position  and  value  of  the  load  in  the  format  position, 
value . 

For  Each  Boundary  Condition 

7.  Enter  the  numerical  code  of  the  known  boundary  value. 

8.  Enter  the  value  of  the  known  boundary  condition. 

Following  the  problem  input,  the  program  prints  on  the  screen  the 
title  of  the  major  subroutines  it  enters.  Within  the  matrix  inversion 
subroutine  an  accuracy  test  is  performed.  (For  an  explanation  of  the 
individual  subroutines  see  the  METHODOLOGY  section  of  this  appendix.) 


I.  PROGRAM  INPUT  (Continued) 

B.  Response  Input 

1.  Enter  the  number  of  response  points  where  shear,  moment, 
displacement,  and  rotation  are  to  be  calculated. 

For  Each  Response  Point 

2.  Enter  the  position  of  the  response  point 

The  response  subroutine  prints  the  number  of  the  response  point 
where  the  solution  is  being  calculated.  Following  the  last  response 
point  the  printing  of  the  output  begins. 


II.  PROGRAM  OUTPUT 


The  user  input  and  response  output  are  always  printed.  The  user 
has  the  option  to  obtain  a  "DEBUG  PRINTOUT."  The  "DEBUG  PRINTOUT" 
consists  of  a  list  of  the  concentrated  loads  that  approximate  the  con¬ 
tinuous  loads  and  the  boundary  loads  (boundary  moments  are  not  included) . 
The  flexibility  matrix  and  its  inverse  are  also  printed.  Appendix  C 
contains  a  few  printouts.  The  last  problem  includes  an  annotated  copy 
of  the  user  interaction. 
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BIM1D  EXAMPLE  PROBLEMS,  BASIC  VERSION 
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1  10  4:  ON  N  GOSUB  1410,1400,1420,1400:  GOSUB  1430: 


1850  FOR  N  »  0  TO  NK  -  1 

1860  PRINT  "RESPONSE  POINT  "N  *  1"  POSITION  INPUT  ”  ";XA(N) 
1870  IF  XA(N) <  0  OR  XA(N)  >  L  GOTO  1860 
1880  NEXT 
1890  RETURN 
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BEAM  ON  AN  ELASTIC  FOUNDATION 
INDIRECT  BOUNDARY  INTEGRAL  METHOD 
THEORY!  DR.  TED  SHUGAR 
PROGRAM:  JAMES  V.  COX 


USER  INPUT 

SPRING 

MODULUS  OF 

MOMENT  OF 

CONSTANT 

ELASTICITY 

INERTIA 

SPAN 

1000. 00 

10000000.00 

25.00 

SO .  00 

CONCENTRATED 
LOAD  NO. 

POSITION 

VALUE 

1 

20.  00 

1 00000 . 00 

2 

60. 00 

- 1 00000 . 00 

LEFT  END 

RIGHT 

END 

MOMENT 

=  0.00 

MOMENT 

S2 

0.  00 

SHEAR 

=  0 . 00 

SHEAR 

‘ 

0.  00 

OUTPUT 

RESPONSE 
POINT  NO. 

POSITION 

SHEAR 

MOMENT 

DISPLACEMENT 

ROTATION 

1 

.  01 

0.00 

0 . 00 

3.506 

-.  070 

■7 

4b. 

10.00 

31536.93 

163395.97 

2.803 

-.072 

3 

20.01 

-44185.37 

606177.57 

2.024 

-.087 

4 

30.00 

-28616.21 

250467. 73 

1 . 064 

-.  103 

5 

40.00 

-23254. 13 

0 . 00 

0. 000 

-.  108 

6 

50.00 

-28616.21 

-250467.73 

-1.064 

-.  103 

7 

60.01 

55782.23 

-606084. 80 

-2.026 

-.087 

8 

70.00 

31536.83 

-163395.97 

-2.803 

-.072 

9 

79.99 

0.00 

0.00 

-3.506 

-.070 

C-ll 


CONCENTRATED 

LOAD  NO. 

POSITION 

VALUE 

o 

4 

0 . 00 

80.00 

145799.80 

-145799.80 

G  MATRIX 

0 

1 

0 

7.902 

- .  500 

-.878 

-.033 

1 

- .  500 

.  016 

-.  033 

0.  000 

2 

-.878 

•  033 

7.  902 

.  500 

3 

•  033 

0 . 000 

.  500 

.016 

GI  MATRIX 

0 

1 

2 

0 

128 

-4.086 

9E-03 

-.639 

1 

■4 . 037 

-64.853 

.  1 10 

-13. 050 

2 

9E-03 

.  639 

-.  128 

4.086 

3 

110 

-13.050 

4.037 

-64.853 

BEAM  ora  AN  ELASTIC  FOUNDATION 
INDIRECT  BOUNDARY  INTEGRAL  METHOD 
THEORY:  DR.  TED  SHUGAR 
PROGRAM:  JAMES  V.  COX 


USER  INPUT 

SPRING  MODULUS  OF  MOMENT  OF 

CONSTANT  ELASTICITY  INERTIA  SPAN 

1000.00  10000000.00  25.00  40.00 

CONCENTRATED 

LOAD  NO.  POSITION  VALUE 


1  20 . 00  1 00000 . 00 

LEFT  END  RIGHT  END 


MOMENT  =  0.00  DISPLACEMENT  =  0.00 

SHEAR  =  0.00  MOMENT  -  0.00 


OUTPUT 


RESPONSE 
POINT  NO. 

POSITION 

SHEAR 

MOMENT 

DISPLACEMENT 

ROTATION 

1 

0.00 

0.  00 

0.00 

3.505 

- .  070 

2 

10.00 

31541.23 

163486.32 

2.802 

-.072 

3 

20.00 

-44196.82 

606448.52 

2.024 

-.087 

4 

30.00 

-28625.78 

250496.21 

1 . 064 

103 

5 

40.  00 

—23268. 66 

0.00 

0.000 

-.  108 
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CONCENTRATED 


LOAD  NO. 

P 

OSITION 

VALUE 

V 

0.  00 

145741.48 

40.00 

-77289.86 

G  MATRIX 

0 

1 

n 

0 

7.904 

-.500 

- 1 . 456 

1 

- .  500 

.  0 1 6 

.  043 

2 

0 . 000 

0 . 000 

0.000 

v>  **" 

■1.456 

-.043 

7.904 

GI  MATRIX 

0 

1 

2 

0 

-.  133 

-4.724 

-32267. 163 

1 

4.  146 

-77.903 

-784819.603 

2 

-.022 

.347 

61323.234 

.ji 

-.  408 

-25.876 

-1 130194.410 

.  043 
6E-03 
0 . 000 
-  500 


.  06 
1.22 
-2E-0 
2.325 
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BEAM  ON  AN  ELASTIC  FOUNDATION 
INDIRECT  BOUNDARY  INTEGRAL  METHOD 
THEORY:  DR.  TED  SHUGAR 
PROGRAM:  JAMES  V.  COX 


CONCENTRATED 


NO. 

POSITION 

VALUE 

1 

2 . 00 

76.  00 

9 

6 . 00 

68 . 00 

.  ji 

1 0 . 00 

60 . 00 

4 

1 4 . 00 

52.00 

5 

18.00 

44 . 00 

6 

22.00 

36. 00 

7 

26.00 

28.00 

8 

30 . 00 

20 . 00 

9 

34.00 

1 2 . 00 

10 

38 . 00 

4 . 00 

11 

42. 00 

-4 . 00 

12. 

46.00 

-12. 00 

13 

50 .  00 

—20. 00 

14 

54.00 

-28 00 

15 

58 . 00 

-36. 00 

16 

62 . 00 

-  44 ..  00 

17 

66 . 00 

-52 . 00 

18 

70 . 00 

-60 . 00 

19 

74 . 00 

-68.00 

20 

78.00 

-76 00 

21 

0 .  00 

880. 14 

22 

80 . 00 

-880. 14 

G  MATRIX 

<; 

;> 

1 

0 

7 . 902 

- .  500 

l 

-.500 

.  016 

2 

-.878 

•  033 

0* 

.  033 

0.  000 

-  .  878 


7 . 902 
.  500 


033 
0 . 000 
.  500 
.016 


GI  MATRIX 

O 


0 

1 


-.128 
-4.037 
9E-03 
1 10 


1 

-4.086 
-64.853 
.639 
-13. 050 


9E-03  -.639 

.110  -13.050 

-.128  4.086 

4.037  —64.853 
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BEAM  ON  AN  ELASTIC  FOUNDATION 
INDIRECT  BOUNDARY  INTEGRAL  METHOD 
THEORY:  DR.  TED  SHIJGAR 
PROGRAM:  JAMES  V.  COX 

USER  INPUT 


SPRING  MODULUS  OF  MOMENT  OF 

CONSTANT  ELASTICITY  INERTIA  SPAN 

1 OOO . 00  1 0000000 . 00  25 . 00  40 . 00 

CONTINUOUS  LEFT  END  RIGHT  END  NO.  OF  LOADS 

LOAD  NO.  POSITION  VALUE  POSITION  VALUE  TO  APPROXIMATE 


1  0.00  20.00  40.00  0.00  10 

LEFT  END  RIGHT  END 


MOMENT  =  0.00  DISPLACEMENT  =  0.00 

SHEAR  -  0 . 00  MOMENT  =  0 . 00 


OUTPUT 


RESPONSE 

POINT  NO. 

POSITION 

SHEAR 

MOMENT 

DISPLACEMENT 

ROTATION 

1 

0.  00 

0 . 00 

0.  00 

.  020 

0  -  000 

10.  00 

-29.43 

33.52 

.015 

0 . 000 

.j. 

20 . 00 

-.74 

4.01 

.  010 

0.  000 

4 

30.  00 

-9.86 

12.43 

5E-03 

0 . 000 

5 

40.  00 

-.94 

0 . 00 

0 . 000 

0 . 000 

CONCENTRATED 


LOAD  NO. 

POSITION 

VALUE 

1 

2 . 00 

76.00 

2 

6 . 00 

68 . 00 

3 

10.  00 

60.  00 

4 

1 4 . 00 

52.00 

5 

18.  00 

44 . 00 

6 

22.00 

36 . 00 

7 

26.  00 

28.00 

8 

30.00 

20.00 

9 

34.00 

12.00 

10 

38.00 

4.00 

11 

0 . 00 

879.82 

12 

40 . 00 

-250.66 

G  MATRIX 

0 

1 

2 

0 

7.904 

-.500 

-1.456 

.  043 

1 

-.500 

.016 

.  043 

6E— 03 

2 

0.000 

0 . 000 

0.000 

0.000 

.J.i  — 

-1.456 

-.043 

7.904 

.  500 

GI  MATRIX 

0 

1 

2 

3 

0 

-.  138 

-4.724 

-32267. 163 

.065 

1 

•4.  146 

-77.903 

-784819.603 

1.226 

2 

-.022 

.347 

61323.234 

-2E-03 

3 

-.408 

-25.876 

-1130194.410 

2.325 

1 . 


Annotated  User  Interaction 


3  RUN 

BOUNDARY  INTEGRAL  METHOD 
SPRING  CONSTANT  =  1000 
E  =  10E6 
I  =  25. 

L.  =■  1 00 

NUMBER  OF  CONTINUOUS  LOADS  =  ± 

CONTINUOUS  LOAD  NUMBER  1 
POSITION  AND  VALUE  OF  LOAD  ON 
LEFT  END  =  20.2000 
POSITION  AND  VALUE  OF  LOAD  ON 
RIGHT  END  =  60,500 

NUMBER  OF  CONCENTRATED  L  "'ADS  TO  APPROX¬ 
IMATE  THE  CONTINUOUS  LOAl- 

NUMBER  OF  CONCENTRATED  LOADS  -  JL 

CONCENTRATED  LOAD  NUMBER  1 

POST I ON  AND  VALUE  OF  LOAD  =  80. 10E5 

BOUNDARY  CONDITION  INPUT 

ENTER  APPROPRIATE  NUMERIC  CODE  FOR  KNOWN 
BOUNDARY  CONDITION 

1.  DISPLACEMENT 

2.  ROTATION 

3.  MOMENT 

4.  SHEAR 


•RESPONSE  INPUT 

NUMBER  OF  RESPONSE  POINTS 
RESPONSE  POINT  1  POSITION 
C> 

RESPONSE  POINT  2  POSITION 
10 

RESPONSE  POINT  3  POSITION 

20 

RESPONSE  POINT  4  POSITION 

30 

RESPONSE  POINT  5  POSITION 
40 

RESPONSE  POINT  6  POSITION 
50 

RESPONSE  POINT  7  POSITION 
60 

RESPONSE  POINT  B  POSITION 
70 

RESPONSE  POINT  9  POSIT I UN 
80 

RESPONSE  POINT  10  POSITION 
90 

RESPONSE  POINT  11  POSITION 
100 

•  RESPONSE 


LEFT  BOUNDARY 

BOUNDARY  CONDITION  1,  CODE  =  ± 
DISPLACEMENT  VALUE  -  0 


CODE 


BOUNDARY  CONDITION  2, 
MOMENT  VALUE  =  0. 

RIGHT  BOUNDARY 


BOUNDARY  CONDITION  3,  CODE  «  _2 
ROTATION  VALUE  =  0. 

BOUNDARY  CONDITION  4,  CODE  =  ± 
SHEAR  VALUE  = 

•GENERATING  CONCENTRATED  LOADS 
•MATRIX  DEVELOPEMENT 
G  MATRIX 
HP  MATRIX 
•MATRIX  INVERSION 
INVERSION  TEST  =  4.17263406 
•SOLVE  FOR  ACTUAL  BEAM  UNKNOWNS 


RESPONSE 

RESPONSE 

RESPONSE 

RESPONSE 

RESPONSE 

RESPONSE 

RESPONSE 

RESPONSE 

RESPONSE 

RESPONSE 

RESPONSE 


PT  1 
PT  2 
PT  3 
PT  4 
PT  5 
PT  6 
PT  7 
PT  8 
PT  9 
PT  10 
PT  11 


Underfilled  quantities  represent  user  input 
•  ■  indicates  die  printout  of  a  subroutine 


BEAM  ON  AN  ELASTIC  FOUNDATION 
INDIRECT  BOUNDARY  INTEGRAL.  METHOD 
THEORY:  DR.  TED  SHIJGAR 
PROGRAM:  TAMES  V.  COX 


USER  INPUT 


, — 

4777777777 


SPRING 

MODULUS  OF 

MOMENT  OF 

CONSTANT 

ELASTICITY 

INERTIA 

SPAN 

1000. 00 

10000000. 00 

25.  00 

1 00 . 00 

CONTINUOUS 

LEFT  END 

RIGHT  END 

NO.  OF  LOADS 

LOAD  NO. 

POSITION  VALUE 

POSITION 

VALUE 

TO  APPROX  IMA 

1 

20.00  2000.00 

60 . 00 

500. 00 

20 

CONCENTRATED 

LOAD  NO. 

POSITION 

VALUE 

1 

80.  00 

1 0000 . 00 

LEFT  END 

RIGHT  END 

DISPLACEMENT 

=  0 . 00 

ROTATION 

SS 

0 . 00 

MOMENT 

=  0.00 

SHEAR 

= 

0 . 00 

OUTPUT 


RESPONSE 
POINT  NO. 

POSITION 

SHEAR 

1 

.01 

7596.64 

2 

10.00 

8965.02 

T 

20.00 

12914.83 

4 

30.00 

790.72 

5 

40.00 

-6419.81 

6 

50.00 

-9662.24 

7 

60.  00 

-9651.91 

9 

70.00 

-3616.04 

9 

80.01 

-8444.32 

10 

90.00 

-4005.30 

11 

99.99 

0.  00 

MOMENT 

DISPLACEMENT 

ROTATION 

0.00 

0.  000 

.028 

80464.25 

.271 

.026 

187869.61 

.510 

.021 

252020. 07 

.675 

.012 

220317. 12 

.742 

2E-03 

137073. 14 

.722 

-5E-03 

38122.64 

.649 

-9E-03 

-27461.57 

.558 

-9E-03 

-37195.25 

.477 

-7E-03 

-98869.39 

.416 

-4E-03 

-1 18683. 19 

.393 

0.  000 

DO  YOU  WANT  THE  DEBUG  PRINTOUT  <Y/N)?Y 


CONCENTRATED 

LOAD  NO.  1 

POSITION 

VALUE 

*? 

21.00 

3925 . 00 

c» 

23 . 00 

3775.00 

4 

25.  00 

3625. 00 

5 

27 . 00 

3475. 00 

6 

29  „  00 

3325. 00 

7 

3 1 . 00 

3175.00 

8 

33.  00 

3025. 00 

9 

35.  00 

2875.00 

10 

37.00 

2725. 00 

11 

39.  00 

2575.00 

12 

4 1 . 00 

2425.00 

13 

43.  00 

2275 . 00 

14 

45.  00 

2125. 00 

15 

47.00 

1975.00 

16 

49 . 00 

1825.00 

17 

5 1 . 00 

1675.00 

18 

53.  O0 

1 525 . 00 

19 

55 . 00 

1375.00 

20 

57 . 00 

1225.00 

21 

59.00 

1075.00 

22 

0 . 00 

-21444.46 

23 

1 00 . 00 

12439.27 

G  MATRIX 

0 

1 

2 

0* 

0  0.000 

0.000 

0 . 000 

0.000 

1  7/901 

-.  500 

-.328 

-.021 

2  0.000 

0.000 

0.000 

0.  000 

3  .021 

-IE-03 

.  500 

.016 

61  MATRIX 

0 

1 

2 

3 

0  63152.345 

0.000 

41453.797 

.086 

1  998138.796 

-1.998 

1337181.850 

.053 

2  -2674.363 

0.000 

998913.625 

1.996 

3  43117.033 

-.083- 

31576172.600 

.025 

Appendix  D 
BIM2D  LISTINGS 


o  o  o  o  n  n  no  ononno  o  onn 


PROGRAM  BIM2D ( INPUT , OUTPUT , TAPE5= INPUT , TAPE6=0UTPUT , TAPE7 ) 
MAIN  PROGRAM 

DIMENSION  X(I52) ,Y(152) ,BC(304) ,KBC(152) ,P(304) ,LN0DE(5) 
REAL  MU, K( 304, 304), LAMDA 
COMMON  /CEES/  CO, Cl ,C2,C3,C4 
MK  =  304 

INPUT  AND  INPUT  PRINT  ROUTINES 

CALL  INFO (X , Y , KBC , BC , NOE , E , MU , KPROB , LNODE , MK ) 

CALL  INPRNT (X , Y ,KBC , BC , NOE , E ,MU , KPROB , LNODE) 

PI  =  2.*ACOS(0. ) 

CALCULATIONS  WHICH  ARE  FUNCTIONS  OF  THE  MATERIAL  ONLY 

G  =  E/2./ (l.+MU) 

PLANE  STRAIN 

LAMDA  ■=  E*MU/(1.+MU)/(1.-2.*MU) 

PLANE  STRESS 

IF (KPROB. EQ.O)  LAMDA  =  2 . *LAMDA*G/ (LAMDA+2 . *G) 

CO  =  1.  +  2 . * (LAMDA+G) /G 
C4  =  PI* (LAMDA  +  2.*G) 

Cl  -  G/2./C4 

C2  =  (LAMDA+G) /C4 

C3  =  (LAMDA+3 . *G) /4 . /C4/G 

C4  =  (LAMDA+G) /4./G/C4 

INFLUENCE  MATRIX  DEVELOPMENT  ROUTINE 

CALL  KMAKER(NOE,LAMDA,G, PI, X,Y,K, KBC, MK, LNODE) 


SOLVE  FOR  ARTIFICIAL  BOUNDARY  LOADS 
N  -  2*NOE 

CALL  AXEQB(K,P,BC,N,MK) 


CALCULATE  RESPONSE  AT  SPECIFIED  POINTS 
CALL  RESPON (NOE , LAMDA ,G,PI,X,Y,P, LNODE) 


SUBROUTINE  MANUAL 


C 

C 

C 

C 

C 

C 

C 

C 

C 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 


VARIABLE  LIST 


A:  A  VALUE  FOR  NONSINGULAR  INTEGRATIONS. 

B:  "  " 

BC(I),BC(I1):  BOUNDARY  TRACTIONS  OR  DISPLACEMENTS  FOR  ELEMENT  (I+I)/2 


c 

CTHETM:  COS(THETAM) . 

c 

CTHETN :  COS (THETAN) . 

c 

CO: 

A  CONSTANT  FOR  A  GIVEN  MATERIAL. 

c 

Cl: 

ti  it 

c 

C2: 

ii  n 

c 

C3: 

it  ti 

c 

C4: 

««  ii 

c 

CSX: 

.EQ.O:  CALCULATE  SX  RESPONSE  AT  A  GIVEN 

POINT 

c 

CSY : 

"  "  SY 

c 

CT: 

II  It  IJ» 

c 

CU: 

..  ..  U 

c 

CV: 

II  II  v 

c 

DL: 

ELEMENT  LENGTH. 

c 

DX: 

FOR  A  GIVEN  ELEMENT  -  X(N)  -  X(N+1). 

c 

DY: 

FOR  A  GIVEN  ELEMENT  -  Y(N)  -  Y(N+1). 

c 

E: 

MODULUS  OF  ELASTICITY. 

c 

EC: 

SQUARE  OF  DISTANCE  FROM  CENTER  OF  M  TO  NODE  1  l 

c 

EE: 

X  DISTANCE  FROM  CENTER  OF  M  TO  NODE  1 

OF 

N. 

c 

EG: 

Y  DISTANCE  FROM  CENTER  OF  M  TO  NODE  1 

OF 

N. 

c 

EPSIL:  A  TOLERANCE  USED  TO  DECIDE  BETWEEN 

INTEGRAL 

c 

G: 

SHEAR  MODULUS  OF  ELASTICITY. 

OF  N. 


GAM1 : 


INTEGRATIONS  OVER  THE  LENGTH  OF  ELEMENT  N. 

II  It 


GAM10:  "  " 

I:  (R)  FIRST  SUBSCRIPT. 

II:  1+1 

J:  (K)  SECOND  SUBSCRIPT. 

Jl:  J+l 

JB:  FLAG  TO  INDICATE  A  BOUNDARY  RESPONSE  POINT. 

K(I,J):  INFLUENCE  MATRIX  TO  OBTAIN  ARTIFICIAL  BOUNDARY  LOADS. 

KBC(M):  KNOWN  BOUNDARY  QUANTITIES  .EQ.O:  STRESSES  AND  .NE.O:  DISPLACEMENTS. 
KPROB:  .EQ.O:  PLANE  STRESS  AND  .NE.O:  PLANE  STRAIN. 

LAMDA:  LAME'  CONSTANT. 

LNODE(NBOUND) :  LAST  NODE  OF  BOUNDARY  NBOUND. 

M:  OUTER  LOOP  ELEMENT  NUMBER  TO  DEVELOP  2  (K)  ROWS. 

MBOUND:  A  BOUNDARY  NUMBER  COUNTER. 

Ml:  M+l 

MR:  MAXIMUM  ORDER  OF  (K) ,  NOE  *  2 

MU:  POISSON'S  RATIO. 
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N:  INNER  LOOP  ELEMENT  NUMBER  TO  DEVELOP  2  (K)  COLUMNS. 

NBOUND:  A  BOUNDARY  NUMBER  COUNTER. 

Nl:  N+l 

NOE:  NUMBER  OF  BOUNDARY  ELEMENTS. 

P(I) ,P(I1) :  ARTIFICIAL  BOUNDARY  LOADS  AT  ELEMENT  (I+l)/2 
PI:  3.14159... 

THETAM:  ANGLE  OF  ELEMENT  M. 

THETAN:  ANGLE  OF  ELEMENT  N. 

STHETM :  SIN (THETAM). 

STHETN:  SIN (THETAN) . 

SXN:  SX  INFLUENCE  COEF  DUE  TO  NORMAL  BOUNDARY  STRESS. 

SXT:  SX  INFLUENCE  COEF  DUE  TO  TANGENTIAL  BOUNDARY  STRESS. 

SXR:  SX  STRESS  AT  A  GIVEN  RESPONSE  POINT. 

SYN:  SY  INFLUENCE  COEF  DUE  TO  NORMAL  BOUNDARY  STRESS. 

SYT:  SY  INFLUENCE  COEF  DUE  TO  TANGENTIAL  BOUNDARY  STRESS. 

SYR:  SY  STRESS  AT  A  GIVEN  RESPONSE  POINT. 

TN:  T  INFLUENCE  COEF  DUE  TO  NORMAL  BOUNDARY  STRESS. 

TT:  T  INFLUENCE  COEF  DUE  TO  TANGENTIAL  BOUNDARY  STRESS. 

TR:  T  STRESS  AT  A  GIVEN  RESPONSE  POINT. 

UN:  U  INFLUENCE  COEF  DUE  TO  NORMAL  BOUNDARY  STRESS. 

UP:  UPPER  LIMIT  ON  ELEMENT  INTEGRATION. 

UT:  U  INFLUENCE  COEF  DUE  TO  TANGENTIAL  BOUNDARY  STRESS. 

UR:  U  DISPLACEMENT  AT  A  GIVEN  RESPONSE  POINT. 

VN:  V  INFLUENCE  COEF  DUE  TO  NORMAL  BOUNDARY  STRESS. 

VT:  V  INFLUENCE  COEF  DUE  TO  TANGENTIAL  BOUNDARY  STRESS. 

VR:  V  DISPLACEMENT  AT  A  GIVEN  RESPONSE  POINT. 

X(N) :  X  COORDINATE  OF  NODE  N,  THE  FIRST  NODE  OF  ELEMENT  N. 

XQM:  X  COORDINATE  OF  THE  CENTER  OF  ELEMENT  M. 

XR:  X  COORDINATE  OF  THE  RESPONSE  POINT. 

Y (N) :  Y  COORDINATE  OF  NODE  N,  THE  FIRST  NODE  OF  ELEMENT  1J. 

YQM:  Y  COORDINATE  OF  THE  CENTER  OF  ELEMENT  M. 

YR:  Y  COORDINATE  OF  THE  RESPONSE  POINT. 


INPUT  FILE 

THE  FOLLOWING  DESCRIPTION  OF  THE  INPUT  FILE  IS  WRITTEN  AS  THOUGH 
A  DECK  OF  CARDS  IS  BEING  USED  AS  THE  INPUT  MEDIUM.  ONE  MUST 
REFER  TO  THE  VARIABLE  LIST  TO  USE  THIS  BREIF  DESCRIPTION. 


CARD/SET  DATA  DESCRIPTION 


FORMAT 


KPROB  (PROBLEM  TYPE:  PL  STRESS/PL  STRAIN)  (15) 


E,MU  (MATERIAL  PROPERTIES) 


(E10.F10) 


LNODE(l) , . . . ,LNODE(5)  (LAST  NODE  NUMBERS)  (515) 


BOUNDARY  DESCRIPTION  CARDS 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


SET  4 


X(1),Y(1),KBC(1),BC(1),BC(2) 


(2F10.I5.2F10) 


X(NOE) ,Y (NOE) ,KBC(NOE) ,BC(2*N0E-1) ,BC(2*N0E) 


SET  5 


RESPONSE  SPECIFICATION  CARDS 

IF  THE  USER  INDICATES  THAT  THE  RESPONSE  IS  BEING 
CALCULATED  ON  A  BOUNDARY  ELEMENT  (JB.NE.O)  THEN 
THE  ELEMENT  NUMBER  (I)  MUST  BE  ENTERED,  AND  BOTH 
XR  AND  YR  ARE  SET  TO  THE  MID-ELEMENT  COORDINATES. 
THE  PROGRAM  DOES  NOT  LIMIT  THE  NUMBER  OF  RESPONSE 
POINTS. 


JB,I,,XR,YR,CSX,CSY,CT,CU,CV 


(2I5.2F10.5I5) 


JB,I, ,XR,YR,CSX,CSY,CT,CU,CV 


SUBROUTINE  LIST 


THE  FOLLOWING  SUBROUTINE  LIST  INDICATES  THE  LEVELS  OF  HIERARCHY 
WITHIN  THE  PROGRAM,  BUT  DOES  NOT  COMPLETELY  DEFINE  THE  FLOW 
OF  EXECUTION.  SOME  OF  THE  SECOND  ORDER,  AND  LOWER  LEVEL, 
ROUTINES  ARE  NOT  ALWAYS  EXECUTED.  FUNCTIONS  ARE  NOT  LISTED. 


BIM2D:  MAIN  PROGRAM 


INFO:  INPUTS  DATA,  EXCEPT  FOR  RESPONSE  DATA. 
INPRNT :  PRINTS  THE  INPUT  DATA. 


KMAKER:  GENERATES  THE  INFLUENCE  COEFFICIENT  MATRIX  (K) . 

PREK:  CALCULATES  THETAM,  AND  MATRIX  ROW  NUMBERS  I  AND  II. 
PREK:  CALCULATES  THETAN,  AND  MATRIX  COL  NUMBERS  J  AND  Jl. 
PREGAM:  CALCULATES  VALUES  IN  PREPARATION  FOR  NONSINGULAR 
INTEGRATIONS,  UP,  EE,  EG,  A,  B,  AND  EC. 

GAMMAF:  CALCULATES  NONSINGULAR  INTEGRATIONS  WHEN  ELEMENT 
M  BOUNDARY  STRESSES  ARE  KNOWN. 

SXC:  CONTROL  ROUTINE  FOR  SXN  AND  SXT  CALCULATIONS  SENDS 
CONTROL  TO  SINGULAR  OR  NONSINGULAR  ROUTINE. 

SXNS:  CALCULATES  NONSINGULAR  SXN  AND  SXT. 

SXS:  CALCULATES  SINGULAR  SXN  AND  SXT. 

SYC:  SY  CONTROL  ROUTINE. 

SYNS:  NONSINGULAR  CALCULATIONS. 

SINGULAR  CALCULATIONS 


SYS: 


TC:  T  CONTROL  ROUTINE. 
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•.TTV'r»  :»  ".'it 


c 


TNS:  NONSINGULAR  CALCULATIONS. 

TS:  SINGULAR  CALCULATIONS. 

PREK:  CALCULATES  THETAN,  AND  MATRIX  COL  NUMBERS  J  AND  Jl. 
PREGAM:  CALCULATES  VALUES  IN  PREPARATION  FOR  NONSINGULAR 
INTEGRATIONS,  UP,  EE,  EG,  A,  B,  AND  EC. 

GAMMAU:  CALCULATES  NONSINGULAR  INTEGRATIONS  WHEN  ELEMENT 
M  BOUNDARY  DISPLACEMENTS  ARE  KNOWN. 

UC:  U  CONTROL  ROUTINE. 

UNS:  NONSINGULAR  CALCULATIONS. 

US:  SINGULAR  CALCULATIONS. 

VC:  V  CONTROL  ROUTINE. 

VNS:  NONSINGULAR  CALCULATIONS. 

VS:  SINGULAR  CALCULATIONS. 

AXEQB:  MATRIX  SOLUTION  ROUTINE  TO  SOLVE  FOR  ARTIFICIAL 
BOUNDARY  LOADS. 

FACTOR 

SUBST 

RESPON:  READS  INPUT  DATA  FOR  RESPONSE  CALCULATIONS, 
CALCULATES  THE  USER  SPECIFIED  RESPONSES, 

AND  PRINTS  THE  CALCULATED  RESPONSES. 

THE  ROUTINES  CALLED  FROM  RESPON  ARE  NOT  LISTED.  THEY  ARE 
MUCH  THE  SAME  AS  THOSE  CALLED  FROM  KMAKER. 


RETURN 

END 


V 
k 
s 

V 

4 


t 
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SUBROUTINE  INFO (X , Y , KBC , BC , NOE , E , MU , KPROB , LNODE , MK) 

C 

C  THIS  ROUTINE  READS  ALL  THE  INPUT  DATA  EXCEPT  FOR  THE  RESPONSE 
C  POINT  DATA.  A  TRAILER  CARD  IS  USED  TO  INDICATE  THE  END  OF  THE 
C  ELEMENT  INPUT,  AND  THE  NUMBER  OF  ELEMENTS  IS  COUNTED  TO 
C  PREVENT  TOO  MANY  ELEMENTS. 

C 

DIMENSION  X(l) ,Y(i) ,KBC(1) ,BC(1) ,LNODE(l) 

REAL  MU 
C 

C  ENTER  PROBLEM  TYPE,  PLANE  STRAIN  OR  PLANE  STRESS. 

C 

READ (5, 1000)  KPROB 
1000  FORMAT(I5) 

C 

C  ENTER  MATERIAL  PROPERTIES 
C 

READ(5 , 1010)  E,MU 
1010  FORMAT(E10.3,F10.3) 

C 

C  ENTER  LAST  NODE  NUMBER  FOR  EACH  BOUNDARY 
C 

READ (5 , 1050) (LNODE (I) , 1=1,5) 

1050  FORMAT(5I5) 

C 

C  LOOP  FOR  ELEMENT  INPUT 
C 

C  INITIALIZE 

NOE  =  0 
N  -  0 
C 

10  N  -  N  +  1 

I  -  2*N  -  1 

II  =  I  +  1 

C  CHECK  FOR  EXCESSIVE  ELEMENT  INPUT 
IF(N.GT.MK/2)  GO  TO  20 

C 

C  ELEMENT  INPUT 

READ(5,1020)X(N) ,Y(N) ,KBC(N) ,BC(I) ,BC(I1) 

1020  FORMAT(2F10.3,I10,2F10. 3) 

C 

C  CHECK  FOR  TRAILER 

IF(X(N) .GT.7.777E+6)  RETURN 
NOE  =  N 
GO  TO  10 
C 
C 

C  N  .GT.  MK/2 
C 

20  READ (5, 1030) TR 


2/2 


^EiD-A133  142 
UNCLASSIFIED 


AN  INVESTIGATION  OF  THE  INDIRECT  BOUNDARY  ELEMENT 
METHOD  IN  ONE-  AND  TWO-  <U)  NAVAL  CIVIL  ENGINEERING 
LAB  PORT  HUENEME  Cfl  T  A  SHUGAR  ET  AL.  MAV  83 
NCEL-TN-1664  F/G  12/1  NL 


1030  FORMAT (F10. 3) 

IF(TR.GT.7.777E+6)  RETURN 
N-MK/2 

WRITE (6, 1040)  N 

1040  FORMAT( 1H0 , 38HNUMBER  OF  ELEMENTS  EXCEEDS  MAXIMUM  OF  ,12, 1H.) 
STOP 
END 


SUBROUTINE  INPRNT (X , Y , KBC , BC , NOE , E ,MU , KPROB , LNODE ) 

C 

C  THIS  ROUTINE  PRINTS  THE  INPUT  DATA,  NOT  INCLUDING  RESPONSE  INPUT 
C  IN  TABULAR  FORM. 

C 

DIMENSION  X(l) ,Y(1) ,KBC(1) ,BC(1) ,LN0DE(1) 

REAL  MU 
C 

C  PRINT  PROGRAM  TITLE 
WRITE(6,1000) 

1000  FORMAT(lHl) 

1010  FORMAT ( IX, T22 ,36H************************************) 

1020  FORMAT(lX,T22 , 1H*) 

1025  FORMAT ( IX, T22 , 1H* ,T57 , 1H*) 

1030  FORMAT (1H+,T5 7 , 1H*) 

WRITE(6,1010) 

WRITE(6,1025) 

WRITE(6,1025) 

WRITE(6, 1020) 

WRITE(6, 1040) 

1040  FORMAT(lH+,T37 , 5HBIM2D) 

WRITE(6, 1030) 

WRITE(6, 1025) 

WRITE (6, 1020) 

WRITE (6, 1050) 

1050  F0RMAT(1H+,T27 ,27HA  BOUNDARY  INTEGRAL  PROGRAM) 

WRITE(6, 1030) 

WRITE (6, 1020) 

WRITE(6,1060) 

1060  FORMAT ( 1H+ ,T27 , 29HFOR  2D  ELASTOSTATICS  PROBLEMS) 

WRITE(6, 1030) 

WRITE(6, 1025) 

WRITE(6,1020) 

WRITE (6, 1070) 

1070  FORMAT (1H+.T32, 17HDEVELOPED  AT  NCEL) 

WRITE(6, 1030) 

WRITE(6,1025) 

WRITE (6, 1025) 

WRITE(6,1010) 

WRITE(6, 1080) 

1080  FORMAT (//////) 


WRITE (6, 11 10) 

IF(KPROB.EQ.O)  WRITE(6, 1090) 

IF(KPROB.NE.O)  WRITE(6,1100) 

1110  FORMAT (1X,24HPR0BLEM  TYPE:  PLANE  STR) 

1090  FORMAT (1H+,T26,3HESS) 

1100  FORMAT ( 1 H+ , T26 , 3HAIN ) 

C 

C  MATERIAL  PROPERTIES 
WRITE(6, 1120) 

1120  FORMAT (1H0,20»1ATERIAL  PROPERTIES, T30,1HE,T45,2HMU) 
WRITE(6,1130)  E,MU 
1130  FORMAT(1X,T22,E10.3,T39,F10.3) 

C 

C  MULTIPLE  BOUNDARY  LAST  NODES 

WRITE(6, 1200)  (LNODE(I), 1-1,5) 

1200  FORMAT (1 HO, 19HB0UNDARY  LAST  NODES/ 5 15) 

C 

C  ELEMENT  DATA 

WRITE (6, 1140) 

1140  F0RMAT(1H0 , 12HELEMENT  DATA) 

WRITE (6, 1150) 

1150  FORMAT ( IX, T10,3HN0. ,T25, 1HX,T35, 1HY ,T41 ,8HKNOWN  BC,T55,3HBC1 
&T65.3HBC2) 

C 

DO  10  N-l.NOE 
WRITE(6, 1160)N,X(N) ,Y(N) 

1160  FORMAT(1X,I12,F10.3,F10.3) 

IF(KBC(N).EQ.O)  WRITE(6, 1170) 

1170  FORMAT (1H+,T45,1HF) 

IF(KBC(N).NE.O)  WRITE(6 , 1180) 

1180  FORMAT (1H+.T45.1 HU) 

I  -  2*N  -  1 

II  -  I  +  1 

WRITE(6,1190)BC(I),BC(I1) 

1190  FORMAT(1H+,T50,2F10. 3) 

10  CONTINUE 
RETURN 
END 


SUBROUTINE  KMAKER (NOE , LAMDA , G , P I , X , Y , K , KBC ,MK , LNODE) 


MATRIX  GENERATION  ROUTINE 

THIS  ROUTINE  DEVELOPS  (K)  FOR  THE  EQUATION  (K)(P)-(BC). 
WHERE  (P)  IS  THE  THE  COLUMN  MATRIX  OF  ARTIFICIAL  BOUNDARY 
STRESSES  AND  (BC)  IS  THE  COLUMN  MATRIX  OF  KNOWN 
BOUNDARY  CONDITIONS. 


j] 


c 


DIMENSION  X(l),Y(l),KBC(l),LNODE(l) 

REAL  K(MK, 1) ,  LAMDA 

COMMON  ASHOE/  F0R,F1R,F2R,F3R,F4R,F01,F11,F21,F31,F41,F02,F12, 
&F22 ,F32,F42 ,FLG 
COMMON  /CEES/  CO, C1,C2,C3,C4 
C 

C  SET  BOUNDARY  NUMBER 
MBOUND  -  I 
DO  10  M-l.NOE 

IF  (M . GT . LNODE (MBOUND ) )  MBOUND  -  MBOUND  +  1 

CALL  PREK (M, Ml ,1,11, STHETM, CTHETM, DL,X,Y, LNODE, MBOUND) 

XQM  -  (X(M)+X(Ml))/2. 

YQM  -  (Y(M)+Y(Ml))/2. 

N BOUND  -  1 

IF(KBC(M).NE.O)  GO  TO  20 
C  BOUNDARY  STRESSES  KNOWN 

DO  30  N-l.NOE 

IF  (N. GT. LNODE (NBOUND))  NBOUND  -  NBOUND  +  1 
CALL  PREK (N , N1 , J , J 1 , STHETN , CTHETN , DL , X , Y , LNODE , NBOUND ) 
IF(I.EQ.J)  GO  TO  40 
C  FOR  NONSINGULARITY 

CALL  PREGAM(UP,DL,A,B, EC, EE, EG,XQM,YQM,X,Y, CTHETN, 

4STHETN.N) 

CALL  GAMMAF(A,B,EC,UP,GAM1,GAM2,GAM3,GAM4,GAM5,GAM6, EE, CTHETN, 
AEG, STHETN) 

40  CALL  SXC(SXN,SXT, STHETN, CTHETN, GAM1.GAM2, LAMDA, G, I, J) 

CALL  SY C ( SYN , SYT , STHETN .CTHETN , GAM3 , GAM4 , LAMDA , G , I , J ) 

CALL  TC(TN,TT, STHETN, CTHETN, GAM5.GAM6, LAMDA, G, I, J) 

K(I, J)  -  STHETM* STHETM* SXN  +  CTHETM* CTHETM* SYN  - 
&2 . *  S THETM* CTHETM* TN 

K(I, Jl)  -  STHETM* STHETM* SXT  +  CTHETM* CTHETM* SYT  - 
&2 . *STHETM*CTHETM*TT 

K(I1 , J)  -  -CTHETM* STHETM* SXN  +  CTHETM* STHETM* SYN  + 

& (CTHETM* CTHETM  -  STHETM*  STHETM) *TN 
K(I1,J1)  -  -CTHETM* STHETM* SXT  +  CTHETM* STHETM* SYT  + 

& (CTHETM* CTHETM  -  STHETM* STHETM) *TT 
30  CONTINUE 

GO  TO  10 
C 
C 

C  BOUNDARY  DISPLACEMENTS  KNOWN 

20  DO  50  N-l.NOE 

IF  (N . GT . LNODE (NBOUND ) )  NBOUND  -  NBOUND  +  1 
CALL  PRER(N,N1,J,J1, STHETN, CTHETN, DL,X,Y, LNODE, NBOUND) 
IF(I.EQ.J)  GO  TO  60 
C  FOR  NONSINGULARITY 

CALL  PREGAM(UP,DL,A,B,EC,EE,EG,XQM,YQM,X,Y,CTHETN, 

&STHETN.N) 

CALL  GAM1AU(A,B,EC,UP,GAM7 ,GAM8,GAM9,GAM10, EE, CTHETN, EG, STHETN) 
60  CALL  UC(UN, UT, STHETN, CTHETN, GAM7, GAMS ,DL,C3,C4, I, J) 

CALL  VC (VN,VT, STHETN, CTHETN, GAM9,GAM10,DL,C3,C4, I, J) 

K(l, J)  -  UN 
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K(I, Jl)  -  UT 
K(I1,J)  -  VN 
K(I1 , Jl)  =  VT 
50  CONTINUE 
10  CONTINUE 
RETURN 
END 


SUBROUTINE  PREK (N , N 1 , J , J 1 , STHETN , CTHETN , DL , X , Y , LNODE , NBOUND) 

DETERMINES  ELEMENT  ANGLE  AND  SEVERAL  ARRAY  INDICES 

DIMENSION  X(1),Y(1),LN0DE(1) 

N1  -  N+l 

IF  (N.EQ. LNODE (NBOUND). AND. NBOUND. EQ.l)  Nl-1 

IF  (N.EQ. LNODE (NBOUND) .AND. NBOUND. NE.l)  Nl=LNODE(NBOUND-l)+l 

DX  -  X(N1)  -  X(N) 

DY  *  Y(N1)  -  Y(N) 

DL  -  SQRT (DX*DX  +  DY*DY) 

STHETN  -  DY/DL 
CTHETN  -  DX/DL 
J  -  2*N  -  1 
Jl  -  J  +  1 
RETURN 
END 


SUBROUTINE  PREGAM (UP , DL , A , B , EC , EE , EG , XQM , YQM , X , Y , CTHETN , 
ASTHETN.N) 

DETERMINES  SEVERAL  VALUES  IN  PREPARATION  FOR  NONSINGULAR 
INTEGRATIONS 

DIMENSION  X(1),Y(1) 

UP  -  DL 

EE  -  XQM  -  X(N) 

EG  -  YQM  -  Y(N) 

A  -  1. 

B  -  -2.* (CTHETN* EE  +  STHETN*EG) 

EC  -  EE*EE  +  EG*EG 

RETURN 

END 
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SUBROUTINE  GAMMAF(A, B,C,UP,GAM1 , GAM2 , GAM3 , GAM4 , GAM5 , GAM6 , 
&E,F,G,H) 


CALCULATES  NONSINGULAR  INTEGRATIONS  WHEN  BOUNDARY  STRESSES 
ARE  KNOWN. 

COMMON  /SHOE/  F0R,F1R,F2R,F3R,F4R,F01,F11,F21,F31,F41, 

1  F02,F12,F22,F32,F42,FLG 

COMMON  /CEES/  CO, Cl ,C2 ,C3 ,C4 
F01-S01(A,B,C,UP) 

Fll“Sll(A,B,C,UP) 

F02-S02(A,B,C,UP) 

F12”S12(A,B,C,UP) 

F22“S22(A,B,C,UP) 

F32-S32(A,B,C,UP) 

Q1-E*F01-F*F11 

Q2-G*F01-H*F11 

Q3“E*G*G*F02-(2.*E*G*H+G*G*F)*F12+(E*H*H+2.*F*G*H)*F22-F*H*H*F32 

Q4«E*E*G*F02-(2.*E*F*G+E*E*H)*F12+(F*F*G+2.*E*F*H)*F22-F*F*H*F32 

GAM1— C0*C1*Q1+C2*Q3 

GAM2«C1*Q2-C2*Q4 

GAM3— C0*C1*Q2+C2*Q4 

GAM4-C1*Q1-C2*Q3 

GAM5-C1*Q2+C2*Q4 

GAM6-C1*Q1+C2*Q3 

RETURN 

END 


SUBROUTINE  GAMMAU (A, B , C , UP , GAM7 , GAM8 , GAM9 , GAM10, 
&E,F,G,H) 


CALCULATES  NONSINGULAR  INTEGRATIONS  WHEN  BOUNDARY  DISPLACEMENTES 
ARE  KNOWN. 

COMMON  /SHOE/  F0R,F1R,F2R,F3R,F4R,F01,F11,F21,F31,F41, 
l  F02,F12»F22,F32,F42,FLG 
COMMON  /CEES/  CO, Cl ,C2,C3,C4 
F01-S01(A,B,C,UP) 

F11-S11(A,B,C,UP) 

F21*S21 (A,B,C,UP) 

FLG-SLG(A,B,C,UP) 

Q5«G*G*F01-2 . *G*H*F1 1+H*H*F2 1 
Q6"E*E*F01-2.*E*F*F11+F*F*F21 
Q7-E*G*F01-(E*H+F*G)*F11+F*H*F21 
Q8-FLG 

GAM7— C3*Q8  -C4*Q5 

GAM8-C4*Q7 

GAM9— C3*Q8-C4*Q6 

GAM10-C4*Q7 

RETURN 

END 
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FUNCTION  S01(A,B,C,UP) 

COMMON  /SHOE/  F0R,F1R,F2R,F3R,F4R,F01 ,F11 ,F21 ,F31 ,F41 
1  F02.F12 ,F22,F32,F42,FLG 

FUN1 (X) -2 . /R*ATAN ( (2 . *A*X+B) /R) 

FUN2(X)-2./(2.*A*X+B) 

FUN3 (X)-l . /R*ALOG( (2 . *A*X+B-R) / (2 . *A*X+B+R) ) 

BAC-B*B  -  4.*A*C 
R-SQRT(ABS(BAC) ) 

EPSIL  -  1.0E-I0 

IF (BAC . GE . -EPSIL . AND . BAC . LE . EPSIL)  GO  TO  2 
IF(BAC)  1,2,3 

1  S01-FUN1(UP)-FUN1(0.) 

RETURN 

2  S01*FUN2 (UP) -FUN2 (0 , ) 

RETURN 

3  S01*FUN3 (UP) -FUN3 (0 . ) 

RETURN 

END 


FUNCTION  Sll(A,B,C,UP) 

COMMON  /SHOE/  F0R,F1R,F2R,F3R,F4R,F01,F11,F21,F31,F41 
1  F02,F12,F22,F32,F42,FLG 

FUN (X) -1 . / (2 . *A) *ALOG ( A*X*X+B*X+C) 

SI l-FUN(UP) -FUN(0 . ) -B/ (2 . *A)*F01 

RETURN 

END 


FUNCTION  S21(A,B,C,UP) 

COMMON  /SHOE/  F0R,F1R,F2R,F3R,F4R,F01 ,F11 ,F21 ,F31 ,F41 
1  F02,F12,F22,F32,F42,FLG 

FUN(X)-X/A-B/ (2 . *A) *AL0G ( A*X*X+B*X+C) 
S21-FUN(UP)-FUN<0. )+(B*B-2 . *A*C) / (2 . *A*A)*F01 
RETURN 
END 


FUNCTION  S02(A,B,C,UP) 

COMMON  /SHOE/  F0R,F1R,F2R,F3R,F4R,F01 ,F11 ,F21 ,F31 ,F41 
1  F02,F12,F22,F32,F42,FLG 

FUN1 (X)-(2 . *A*X+B) / ( (-BAC) * (A*X*X+B*X+C) ) 

FUN2  (X)— 1 .  /A/A/3 .  /  (B/2 .  /A+X)**3 . 

BAC  -  B*B  -  4 . *A*C 
EPSIL  -  1.0E-10 

IF(BAC.GE. -EPSIL. AND. BAC. LE. EPSIL)  GO  TO  2 
IF(BAC)  1,2,1 

1  S02-FUN1 (UP) -FUN1 (0 . ) +2 . *A/ (-BAC) *F01 
return 

2  S02-FUN2 (UP) -FUN2 (0 . ) 

RETURN 


FUNCTION  S12(A,B,C,UP) 

COMMON  /SHOE/  F0R,F1R,F2R,F3R,F4R,F01 ,F11 ,F21 ,F31 ,F41 
1  F02,F12,F22,F32,F42,FLG 

FUN1  (X)— (2 .  *C+B*X)  /  (  (-BAC)  *  (A*X*X+B*X+C)  ) 

FUN2  (X)— (B/2 . +3 .  *A*X)  /  (3 .  *  A*X+ 1 . 5*B) 

BAC  -  B*B  -  4 . * A*C 
EPSIL  -  1.0E-10 

IF(BAC.GE.-EPSIL.AND.BAC.LE.EPSIL)  GO  TO  2 
IF(BAC)  1,2,1 

1  S12-FUNl(UP)-FUNl(0.)-B/(-BAC)*F01 
RETURN 

2  S12-FUN2(UP)-FUN2(0.) 

RETURN 

END 


FUNCTION  S22(A,B,C,UP) 

COMMON  /SHOE/  F0R,F1R,F2R,F3R,F4R,F01,F11,F21,F31,F41 
1  F02,F12,F22,F32,F42,FLG 

FUN  (X)  —X/  (A*  ( A*X*X+B*X+C  )  ) 

S22-FUN(UP) -FUN (0 . ) +C/A*F02 

RETURN 

END 


FUNCTION  S32(A,B,C,UP) 

COMMON  /SHOE/  F0R,F1R,F2R,F3R,F4R,F01.F11,F21,F31,F41 
1  F02 , FI 2 , F22 , F32 , F42 , FLG 

FUN(X)— X*X/  (A*  (A*X*X+B*X+C)  ) 
S32-FUN(UP)-FUN(0.)+F11/A+C*F12/A 
RETURN 
END 


FUNCTION  SLG(A,B,C,UP) 

COMON  /SHOE/  F0R,F1R,F2R,F3R,F4R,F01,F11,F21,F31,F41 
1  F02,F12,F22,F32,F42,FLG 

FUN(X)"X*ALOG(SQRT(A*X*X+B*X+C) ) 

SLG-FUN(UP)-A*F21-B* . 5*Fl 1 

RETURN 

END 


SUBROUTINE  SXC (SXN ,SXT , STHETN , CTHETN »GAMI ,GAM2,LAMDA,G, I ,  J) 

C 

C  SX  CONTROL  ROUTINE  TO  DECIDE  BETWEEN  SINGULAR  AND  NONSINGULAR 
C  INFLUENCE  COEFFICIENTS. 

C 

REAL  LAMDA 
IF  (I.EQ.J)  GO  TO  10 
C  NONSINGULAR 

CALL  SXNS (SXN, SXT, STHETN, CTHETN, GAM1 ,GAM2) 

RETURN 
C  SINGULAR 

10  CALL  SXS( SXN, SXT, STHETN, CTHETN, LAMDA, G) 

RETURN 

END 


C 

C 

C 


SUBROUTINE  SXNS (SXN, SXT, STHETN, CTHETN, GAM 1 ,GAM2) 

SIGMA  X  INFLUENCE  COEFS  FOR  NONSINGULARITY  CONDITIONS 


SXN  -  -STHETN*GAM1  +  CTHETN *GAM2 
SXT  -  CTHETN*GAM1  ♦  STHETN*GAM2 
RETURN 
END 


C 

C 

C 


SUBROUTINE  SXS ( SXN , SXT , STHETN , CTHETN , LAMDA , G) 

SIGMA  X  INFLUENCE  COEFS  FOR  SINGULARITY  CONDITIONS 


REAL  LAMDA 

SXN  -  (LAMDA  +  2.*G*STHETN*STHETN)  /2.  /(LAMDA  +  2.*G) 

SXT  -  -STHETN*CTHETN 

RETURN 

END 


C 

C 

C 


SUBROUTINE  SYC ( SYN , SYT , STHETN , CTHETN , GAM3 , GAMA , LAMDA , G , I , J ) 
SY  CONTROL  ROUTINE 


REAL  LAMDA 
IF  (I.EQ.J)  GO  TO  10 
C  NONSINGULAR 

CALL  SYNS ( SYN , SYT , STHETN , CTHETN , GAM3 , GAMA ) 
RETURN 
C  SINGULAR 

10  CALL  SYS(SYN, SYT, STHETN, CTHETN, LAMDA, G) 
RETURN 
END 
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SUBROUTINE  SYNS ( SYN , SYT , STHETN .CTHETN , GAM3 , GAM4 ) 

SIGMA  Y  INFLUENCE  COEFS  FOR  NONSINGULARITY  CONDITIONS 

SYN  -  CTHETN*GAM3  -  STHETN*GAM4 
SYT  -  STHETN*GAM3  +  CTHETN* GAM4 
RETURN 
END 


SUBROUTINE  SYS (SYN , SYT , STHETN , CTHETN , LAMDA , G) 

SIGMA  Y  INFLUENCE  COEFS  FOR  SINGULARITY  CONDITIONS 
REAL  LAMDA 

SYN  -  (LAMDA  +  2 . *G*CTHETN*CTHETN )  /2.  /(LAMDA  +  2.*G) 

SYT  -  STHETN* CTHETN 

RETURN 

END 


SUBROUTINE  TC (TN , TT , STHETN , CTHETN , GAM5 , GAM6 , LAMDA , G , I , J ) 

TAU  CONTROL  ROUTINE 

REAL  LAMDA 
IF  (I.EQ.J)  GO  TO  10 
NONSINGULAR 

CALL  TNS (TN , TT , STHETN , CTHETN ,GAM5 , GAM6 ) 

RETURN 

SINGULAR 

D  CALL  TS(TN,TT, STHETN, CTHETN, LAMDA.G) 

RETURN 

END 


SUBROUTINE  TNS (TN.TT, STHETN, CTHETN, GAM5 ,GAM6) 

TAU  INFLUENCE  COEFS  FOR  NONSINGULARITY  CONDITIONS 

TN  -  STHETN *GAM5  -  CTHETN*GAM6 
TT  -  -CTHETN*GAM5  -  STHETN*GAM6 
RETURN 
END 


UUU  U  U  O  U  Ur-,  u  u  u  u  o  o 


SUBROUT INE  TS (TN , TT , STHETN , CTHETN , LAMDA , G) 

TAU  INFLUENCE  COEFS  FOR  SINGULARITY  CONDITIONS 
REAL  LAMDA 

TN  -  -G  /  (LAMDA  +  2 . *G) * STHETN* CTHETN 
TT  -  -.5  *  ( STHETN* STHETN  -  CTHETN* CTHETN) 
RETURN 
END 


SUBROUTINE  UC (UN , UT , STHETN , CTHETN , GAM7 , GAM8 ,DL,C3,C4,I,J) 

U  CONTROL  ROUTINE 

IF  (I.EQ.J)  GO  TO  10 
NONSINGULAR 

CALL  UNS (UN, UT, STHETN, CTHETN, GAM7.GAM8) 

RETURN 

SINGULAR 

0  CALL  US(UN,UT, STHETN, CTHETN, DL,C3,C4) 

RETURN 

END 


SUBROUTINE  UNS (UN, UT, STHETN, CTHETN, GAM7 ,GAM8) 

U  INFLUENCE  COEFS  FOR  NONSINGULARITY  CONDITIONS 

UN  *  -STHETN*GAM7  +  CTHETN*GAM8 
UT  -  CTHETN*GAM7  +  STHETN*GAM8 
RETURN 
END 


SUBROUTINE  US (UN, UT, STHETN, CTHETN, DL,C3,C4) 

U  INFLUENCE  COEFS  FOR  SINGULARITY  CONDITIONS 

UN  -  DL  *  STHETN  *  (C3*(ALOG(DL/2. )  -1.)  +  C4) 
UT  -  -C3  *  (ALOG(DL/2.)  -1.)  *  DL  *  CTHETN 
RETURN 
END 


ofio  non  >—  n  n  non 


SUBROUTINE  VC ( VN , VT , STHETN , CTHETN , GAM9 , GAMI 0 , DL , C3 , C4 , I , J) 

V  CONTROL  ROUTINE 

IF  (I.EQ.J)  GO  TO  10 
NONSINGULAR 

CALL  VNS ( VN , VT , STHETN , CTHETN , GAM9 , GAMI 0) 

RETURN 

SINGULAR 

0  CALL  VS (VN,VT, STHETN, CTHETN, DL,C3,C4) 

RETURN 

END 


SUBROUTINE  VNS (VN , VT , STHETN , CTHETN , GAM9 , GAMI 0) 

V  INFLUENCE  COEFS  FOR  NONSINGULARITY  CONDITIONS 

VN  -  CTHETN*GAM9  -  STHETN*GAM10 
VT  -  STHETN*GAM9  +  CTHETN*GAM10 
RETURN 
END 


SUBROUTINE  VS (VN , VT , STHETN , CTHETN , DL , C3 , C4 ) 

V  INFLUENCE  COEFS  FOR  SINGULARITY  CONDITIONS 

VN  -  -DL  *  CTHETN  *  (C3*(ALOG(DL/2.)  -1.)  +  C4) 
VT  -  -C3  *  (ALOG(DL/2.)  -1.)  *  DL  *  STHETN 
RETURN 
END 


! 
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SUBROUTINE  AXEQB(A,X,B,N,M) 

DIMENSION  A(M, 1) ,  X(l),  B(l),  IPIVOT(304),  D(304),  W(304,304) 

C 

C  THIS  SUBROUTINE  SOLVES  THE  LINEAR  SYSTEM  AX  =  B 

C  IN  THIS  APPLICATION  THE  ARTIFICIAL  BOUNDARY  STRESSES  ARE  SOLVED  FOR 
C 

CALL  FACTOR  (A,A,IPIVOT,D,N,IFLAG,M) 

IF  (I FLAG  .EQ.  1)  GO  TO  10 
WRITE  (6,1000) 

STOP 

10  CONTINUE 

DO  100  I  =  1,N 
100  X(I)  =  0.0 

CALL  SUBST  (A,B,X,IPIVOT,N,M) 

RETURN 

1000  FORMAT (19H1MATRIX  IS  SINGULAR) 

END 


SUBROUTINE  FACTOR ( A , W , IP I VOT , D , N , I FLAG , M) 
DIMENSION  A(M,l),IPIVOT(l),D(l) 

DIMENSION  W (M , 1 ) 

I FLAG  -  1 

C  INITIALIZE  W,  IPIVOT,  D 
DO  10  1=1, N 
IPIVOT(I)  -  I 
ROWMAX  =  0. 

DO  9  J  -  1 ,N 
W(I,J)  =  A(I , J) 

9  ROWMAX  =  AMAX1 (ROWMAX, ABS(W(I,J))) 

IF  (ROWMAX  .EQ.  0.)  GO  TO  999 

10  D(I)  -  ROWMAX 

C  GAUSS  ELIMINATION  WITH  SCALED  PARTIAL  PIVOTING. 
NM1  -  N-l 

IF  (NM1  .EQ.  0)  RETURN 
DO  20  K  «  1 ,NM1 
J  =  K 

KP1  -  K  +  1 
IP  «  IPIVOT (K) 

COLMAX  »  ABS (W(IP,K) )/D(IP) 

DO  11  I  -  KPl ,N 
IP  =  IPIVOT(I) 

AWIKOV  =  ABS(W(IP,K))/D(IP) 

IF  (AWIKOV  .LE.  COLMAX)  GO  TO  11 
COLMAX  -  AWIKOV 
J  -  I 

11  CONTINUE 

IF  (COLMAX  .EQ.  0.)  GO  TO  999 
C 

IPK  =  IPIVOT (J) 

IPIVOT (J)  -  IPIVOT (K) 

IPIVOT (K)  -  IPK 


DO  20  I  -  KPl.N 
IP  -  IPIVOT(I) 

W(IP,K)  -  W(IP,K)/W(IPK,K) 

RATIO  -  -W(IP.K) 

DO  20  J  -  KPl.N 

20  W(IP,J)  -  RATI0*W(IPK, J)  +  W(IP,J) 
IP  (W(IP,N)  .EQ.  0.)  GO  TO  999 
RETURN 

999  I FLAG  -  2 
RETURN 
END 


SUBROUT IN S  &UBST ( W , B , X , IP I VOT , N ,  M) 
DIMENSION  W(M,1),B(1) ,X(I) ,IPIVOT(l) 
IF  (N.GT.l)  GO  TO  10 
X(l)  -  B(I)/W(1,1) 

RETURN 

10  IP  -  IPIVOT(l) 

X(l)  -  B(IP) 

DO  15  K  -  2,N 
IP  -  IPIVOT(K) 

KM1  -  K-l 
SUM  -  0. 

DO  14  J-  l.KMl 

14  SUM  -  W(IP,J)*X(J)  +  SUM 

15  X(K)  -  B(IP)  -  SUM 

X(N)  -  X(N)/W(IP,N) 

K  -  N 

DO  20  NP1MK  -  2,N 
KP1  -  K 
K  «  K-  1 
IP  -  IPIVOT(K) 

SUM  -  0. 

DO  19  J  -  KPl.N 

19  SUM  -  W(IP,J)*X(J)  +  SUM 

20  X(K)  -  (X(K)  -  SUM) /W(IP,K) 

RETURN 

END 


SUBROUTINE  RESPON  (NOE,LAMDA,G,PI,X,Y,P,LNODE) 

C 

C  THIS  ROUTINE  CALCULATES  THE  RESPONSE  AT  EACH  USER  DEFINED 
C  POINT  OF  INTEREST.  THE  RESULTS  ARE  SENT  TO  DEVICE  6,  WHICH 
C  IS  INITIALLY  OUTPUT. 

C  THE  USER  HAS  THE  OPTION  OF  DEFINING  THE  RESPONSE  POINT 
C  AS  WITHIN  OR  ON  THE  BOUNDARY.  FOR  THAT  GIVEN  POINT  HE  ALSO 
C  HAS  THE  OPTION  OF  WHICH  RESPONSES  TO  CALCULATE. 

C 

COMMON/CEES/CO,CI ,C2,C3,C4 
DIMENSION  X(l),Y(l),P(l),LNODE(l) 

INTEGER  CSX,CSY,CT,CU,CV 
REAL  LAMDA 
C 

C  PRINT  OUTPUT  HEADER 
WRITE (6, 3000) 

3000  FORMAT(  1H1  ,T4 .8HRESPONSE/  IX  ,T4 , 8H - ) 

WRITE(6,3010) 

3010  FORMAT ( 1H0 ,T4 , 7HELEMENT , T 1 6 , 1HX , T26 , 1HY , T36 , 2HSX , T46 , 2HSY , T56 , 1HT 
&,T66,1HU,T76,1HV/) 

C 

C  RESPONSE  POINT  CALCULATION  LOOP 
C 

10  READ(5,3015)JB,I,XR,YR,CSX,CSY,CT,CU,CV 
3015  FORMAT(2I5,2F10. 0,515) 

IF(JB.EQ. 77777)  RETURN 
IF(JB.EQ.O)  GO  TO  20 
C  RESPONSE  AT  A  BOUNDARY  ELEMENT 
WRITE(6,3020)I 
3020  FORMAT ( IX, T5, 15) 

II  -  I  +  1 
C 

C  DETERMINE  LAST  NODE  NUMBER 
DO  15  N-1,5 

IF(LNODE(N).EQ.O)  GO  TO  16 
IF(I.EQ.LNODE(N) .AND.N.EQ. 1)  11-1 

15  IF(I.EQ.LNODE(N) .AND.N.NE. 1)  Il-LNODE(N-l)+l 
C 

16  XR  -  (X(I)  +  X(Il))/2. 

YR  -  (Y(I)  +  Y(Il))/2. 

GO  TO  30 

C  RESPONSE  NOT  AT  A  BOUNDARY  ELEMENT 
20  WRITE(6,3030)XR,YR 
3030  FORMAT(1X,T11,2F10.3) 

I  -  0 

C  INITIALIZE  THE  RESPONSES 
30  SXR  -  0. 

SYR  -  0. 

TR  -  0. 

UR  -  0. 

VR  -  0. 


C  INITIALIZE  BOUNDARY  NUMBER 
NBOUND- 1 

C  FOR  EACH  LOAD  (ASSUMES  BOUNDARY  LOADS  ONLY) 

C 

DO  40  N-I.NOE 

IF (N . GT . LNODE (NBOUND) )  NBOUND=NBOUND+ I 

CALL  PREK(N,N1,J,J1, STHETN , CTHETN , DL , X , Y , LNODE , NBOUND ) 

CALL  PREGAM(UP,DL, A, B, EC, EE, EG, XR,YR,X,Y, CTHETN, STHETN, N) 
IF(CSX.NE.O.AND.CSY.NE.O.AND.CT.NE.O)  GO  TO  50 
C  THERE  ARE  STRESS  CALCULATIONS 

CALL  GAMMAF ( A , B , EC , UP , GAM1 , GAM2 , GAM3 , GAM4 , GAMS , GAM6 , EE , CTHETN , 
&EG, STHETN) 

IF(CSX.NE.O)  GO  TO  60 

CALL  SXC ( SXN , SXT , STHETN , CTHETN , GAM1 , GAM2 , LAMDA , G , I , N) 

SXR  -  SXR  +  SXN*P(J)  +  SXT*P(J1) 

60  IF(CSY.NE.O)  GO  TO  70 

CALL  SYC ( SYN , SYT , STHETN , CTHETN , GAM3 , GAM4 , LAMDA , G , I , N) 

SYR  -  SYR  +  SYN*P(J)  +  SYT*P(J1) 

70  IF(CT.NE.O)  GO  TO  50 

CALL  TC (TN , TT , STHETN , CTHETN , GAM5 , GAM6 , LAMDA , G , I , N) 

TR  »  TR  +  TN*P(J)  +  TT*P(J1) 

50  IF(CU.NE.O.AND.CV.NE.O)  GO  TO  40 
C  THERE  ARE  DISPLACEMENT  CALCULATIONS 

CALL  GAMMAU (A, B , EC , UP , GAM7 ,GAM8 , GAM9 , GAM10, EE , CTHETN , EG , STHETN) 
IF(CU.NE.O)  GO  TO  80 

CALL  UC ( UN , UT , STHETN , CTHETN , GAM7 , GAM8 , LAMDA , G , I , N) 

UR  -  UR  +  UN*P(J)  +  UT*P(J1) 

80  IF(CV.NE.O)  GO  TO  40 

CALL  VC(VN,VT, STHETN, CTHETN, GAM9,GAM10,LAMDA,G, I, N) 

VR  -  VR  +  VN*P(J)  +  VT*P(J1) 

40  CONTINUE 

C  INDIVIDUAL  RESPOPNSE  POINT  PRINT  OUT 
IF(CSX.NE.O)  GO  TO  90 
WRITE(6,3400)  SXR 
3400  FORMAT (IH+,T31 ,E10.3) 

90  IF(CSY.NE.O)  GO  TO  100 
WRITE(6,3500)  SYR 
3500  FORMAT(1H+,T41,E10.3) 

100  IF(CT.NE.O)  GO  TO  110 
WRITE(6,3600)  TR 
3600  FORMAT (1H+.T51 ,E10. 3) 

110  IF(CU.NE.O)  GO  TO  120 
WRITE(6,3700)  UR 
3700  FORMAT(1H+,T61,E10.3) 

120  IF(CV.NE.O)  GO  TO  130 
WRITE(6,3800)  VR 
3800  FORMAT(1H+,T71,E10.3) 

C  MAKE  OUTPUT  FILE  (TAPE7)  FOR  PLOTTING. 

130  IF(JB.NE.O)  GO  TO  10 

WRITE(7,4000)  XR,YR,SXR,SYR,TR,UR,VR 
4000  FORMAT(2F10.4,5E15.6) 

GO  TO  10 
END 
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